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ABSTRACT 


FINITE  ELEMENT  ANALYSIS  OF  LASER-INDUCED  DAMAGE  TO 
MECHANICALLY  LOADED  LAMINATED  COMPOSITES  IN  HELICOPTERS,  by 
LTC  Mohammed  J.K  Alghatam,  BDF,  Bahrain,  113  pages. 

This  thesis  examines  the  lethality  of  laser-directed  energy  weapons  in  causing 
structural  failure  to  mechanically  loaded  laminated  composites  in  helicopters.  The 
analysis  is  based  on  a  three-dimensional  numerical  finite  element  model  of  a 
rectangular,  incipient,  orthotropic  compression  panel  with  its  edges  rotationally 
restrained.  The  panel  is  made  of  24”x6"  graphite/epoxy  laminated  composite  where 
the  laminate  consists  of  fourteen  laminas  of  different  ply  orientation  angles  and  two 
different  thicknesses. 

The  data  for  the  degradation  of  the  material  properties  with  increasing 
temperature,  as  well  as  temperatime  dependence  thermal  response  data,  were 
provided  by  the  Naval  Research  Laboratory,  Washington,  D.C.  Laser  irradiation  beam 
strikes  the  center  lines  of  the  panel  with  a  radius  of  3”  and  applied  power  load  of  1 
kilowatt  per  square  centimeter  causing  an  intense  localized  heating  to  the  already 
mechanically  loaded  panel. 

Failure  of  laminated  composites  is  controversial,  even  under  ambient 
temperature.  However,  a  successful  attempt  has  been  conducted  here  to  predict  the 
mechanical  and  thermal  buckling  and  post-buckling  behavior  of  the  heated  panel. 
This  is  further  complicated  by  the  anisotropic  behavior  of  the  composite  material  and 
the  thermal  effect.  Even  without  buckling,  the  problem  is  thermally  and  structurally 
non-linear.  This  work  can  be  considered  as  an  above  preliminary  numerical  analysis. 
It  provides  insight  (to  be  verified  by  experimentation)  into  using  existing  codes  to 
obtain  approximate  solutions  good  enough  for  the  engineering  design,  with  speed  and 
low  cost. 

Should  this  expertise  and  knowledge  be  available,  the  designer  of  directed 
energy  laser  weapon  can  optimize  designs  of  future  weapons,  or  make  efficient  use  of 
existing  ones.  Likewise,  the  helicopter  structural  designer  can  optimize  the  design  of 
future  panels  and  structures,  or  make  efficient  use  of  existing  ones.  In  addition,  the 
designer  can  tailor  a  composite  material  to  meet  a  particular  structural  requirement 
with  little  waste  of  material  capability.  The  ability  to  tailor  a  composite  material  to 
its  job  is  the  biggest  advantage  that  composites  have  over  metallic  or  plastic 
structures.  The  complete  prediction  gives  the  values  of  static  and  thermal  loads, 
stresses,  strains,  and  various  other  heat  transfer  parameters  throughout  the  domain 
of  interest. 

Finite  element  method  superiority  to  other  numerical  analysis  techniques  is 
evident  and  numerical  results  of  this  work  compare  favorably  with  experimental 
results. 
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CHAPTER  1 


INTRODUCTION 


This  work  concerns  the  current  directed  energy  program 
investigating  lethality  of  lasers  and  its  support.  A  laser  load  of  sufficient 
magnitude  can  cause  structural  failure  for  ballistic  missiles,  fixed-wing  planes, 
and  helicopters.  Laser  loading  can  damage  main  structural  components  or 
vulnerable  internal  components,  which  can  result  in  the  failure  of  the  missile 
or  aircraft  mission. 

Therefore,  it  is  necessary  to  understand  the  response  of  laminated 
composite  thin  shell  structure  to  sudden  and  intense  loading.  The  objective  of 
the  approach  is  to  develop  a  three-dimensional  numerical  finite  element  model 
and  predict  the  response  of  mechanically  loaded  and  intensely  localized  spot 
heated  composite  structures  in  helicopters.  The  reason  for  using  numerical 
methods  is  due  to  the  anticipated  and  already  practically  verified  non-linearity 
of  the  thermal  and  mechanical  responses. 
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It  should  be  emphasized  that  the  subject  of  the  failure  of  laminated 
composites  is  controversial,  even  under  ambient  temperatures.  Since  the 
processes  undor  consideration  have  a  very  important  impact  in  directed  energy 
weapon  environment,  they  should  be  dealt  with  effectively.  This  can  be  done 
as  a  result  of  understanding  the  processes  and  their  nature  and  the  ability  to 
predict  them  quantitatively.  If  this  expertise  and  knowledge  is  available,  the 
designer  of  a  directed  energy  laser  weapon  can  ensure  the  desired  performance 
and  can  choose  an  optimum  design  from  many  alternate  possibilities.  Also  the 
power  of  prediction  enables  the  operation  of  existing  laser  weapons  more 
efficiently.  On  the  other  hand,  the  helicopter  structural  designer  can 
determine  critical  laminated  composite  panels  of  the  structure  needing 
strengthening  or  redundancy  to  survive  laser  strikes  at  a  relatively  low  cost. 

The  prediction  of  behavior  in  a  given  physical  situation  requires  the 
values  of  the  relevant  variables  governing  the  process  of  interest.1  In  this 
thesis  work  on  directed  energy  laser  weapons  against  helicopters,  the  complete 
prediction  gives  us  the  values  of  static  and  thermal  loads,  stresses,  and  strains; 
it  also  provides  heat  transfer  parameters  such  as  temperature  distributions 
throughout  the  domain  of  interest. 

Actual  measurement  gives  the  most  reliable  information  about  any 
physical  process.  Since  the  Naval  Research  Laboratory  (NRL)  has  conducted 


'S.  V.  Patankar,  Numerical  Heat  Transfer  and  Fluid  Flow  (New  York:  Hemisphere  Publishing 
Corporation,  McGraw-Hill  Book  Company,  1980),  p.  12. 
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such  experiments,  they  were  kind  enough  to  forward  all  available  physical 
data. 

As  for  the  numerical  method  we  are  going  to  use  the  theoretical 
prediction,  it  works  out  the  results  of  a  mathematical  model  consisting  of  a  set 
of  non-symmetrical,  and  non-linear  partial  differential  equations.  The 
theoretical  calculations  offer  many  advantage  over  a  corresponding 
experimental  investigation.  Such  advantages  are  apparent  in  low  cost,  speed, 
complete  information,  and  ability  to  simulate  realistic  and  ideal  conditions. 
Still,  there  are  drawbacks  and  limitations,  the  most  imoortant  of  which  is  that 
the  validity  of  the  mathematical  model  limits  the  usefulness  of  the 
computation. 

This  work  involves  a  difficult  problem  of  complex  composites 
geometry  and  sensitive  materials  properties,  as  well  as  non-linearities.  The 
extremely  intense  high  temperatures,  incipient  structures  with  post-buckled 
loadings  and  spot  heating,  make  it  hard  to  solve  numerically.  Up  to  now,  to 

the  best  of  the  author’s  knowledge,  no  one  has  analyzed  the  problem  of  laser 

* 

damage  in  loaded  laminated  composite  panels  (as  configured  here)  using  the 
finite  element  method. 

The  finite  element  method  is  a  general  numerical  method  based  on 
an  approximate  solution  of  a  continuum  problem  that  allows  an  immediate 
extension  to  non-structural  problems.  We  define  continuum  as  a  body  of 
matter  (solid,  liquid,  gas,  or  plasma)  or  simply  a  region  of  space  in  which  a 
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particular  phenomenon  is  occurring.2  Quantities  with  obvious  physical 
meaning  are  chosen  as  the  variable  parameters,  hence  giving  a  direct  physical 
contact  with  the  real  problem  being  examined. 

The  label  “finite  element  method*  first  appeared  in  1954,  v  hen  it 
was  used  by  R.W.  Clough3  in  a  paper  on  plane  elasticity  problems,  but  the 
ideas  of  finite  element  analysis  date  back  much  further. 

The  foundation  of  framework  analysis  was  laid  in  the  period 
1850-1875  by  Maxwell  and  Mohr,  amongst  others.  These  concepts  provide  the 
methodology  of  matrix  structural  analysis  that  underlies  finite  element  theory. 
A  lack  of  quick  solution  methods  for  multiple  algebraic  equations  prevented 
any  true  further  advance  for  eighty  years.  These  limitations  were  relieved 
somewhat  in  1932  when  Hardy  Cross  introduced  the  method  of  Moment 
Distribution.  This  made  it  feasible  to  solve  more  complex  problems  by  orders 
of  magnitude. 

The  introduction  of  the  computer  in  the  1950s  gave  a  means  of 
handling  quickly  the  previously  unsolvable  matrix  algebra.  The  solution  of 
real  problems  dealing  with  complex  continua  by  limiting  their  infinite  degrees 
of  freedom  to  a  finite  number  of  unknowns  became  possible.  Such  a  process 


2Mohammed  J.  Alghatam  “An  Investigation  into  the  Use  of  the  Finite  Element  Method  for  Certain 
Continuum  Problems”  (Unpublished  MSC  thesis,  University  of  Technology,  Loughborough,  Leics, 
England,  UK,  1979),  p.  1. 

3Robert  W.  Clough,  “The  Finite  Element  Method  in  Structural  Mechanics,"  in  Stress  Analysis: 
Recent  Developments  in  Numerical  and  Experimental  Methods,  ed.  by  O.C.  Zienkiewicz  and  G.S. 
Holister  (London:  Wiley  and  Sons  Ltd.,  1954),  pp.  48-79. 
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of  discretization  was  first  successfully  performed  by  the  now  well  known 
method  of  “Finite  Differences.” 

The  mid-fifties  saw  the  arrival  of  another  approach,  that  of  finite 
elements.  Its  somewhat  simple  logic  makes  it  ideally  suited  for  the  computer. 
The  first  elements  evolved  were  by  Turner,  Clough,4  and  Melosh,  triangular, 
rectangular,  and  plate  bending  respectively. 

Although  originally  derived  to  study  stress  in  complex  airframe 
structural  problems,  finite  element  methods  have  expanded  into  such  fields  as 
heat  transfer  and  fluid  flow.  This  extension  process  of  the  finite  element 
technique  (the  variational  approach)  to  non-structural  situations,  was 
documented  by  O.C.  Zienkiewicz  and  Cheung  in  1965.5,8,7 

The  finite  element  method,  as  a  numerical  analysis  technique  for 
obtaining  approximate  solutions  to  a  wide  variety  of  engineering  problems,  has 
extend*4  and  been  applied  to  the  broad  field  of  continuum  mechanics.  The 
finite  element  method  has  received  much  attention  in  engineering  schools  and 
industry  due  to  its  diversity  and  flexibility  as  an  analytical  tool. 


*Ibid.  pp.  48-79. 

80.  C.  Zienkiewicz,  The  Finite  Element  Method  (London:  McGraw-Hill  Book  Company,  1977),  p. 
435. 

®Y.  K.  Cheung  and  M.  F.  Yeo,  Practical  Introduction  to  Finite  Element  AnalvBiB  (London:  Pitman 
Publishing  Co.,  1979),  p.  112. 

70.  C.  Zienkiewicz,  “The  Finite  Element  Method:  From  Intuition  to  Generality.”  Applied  Mechanics 
Review  (London:  McGraw-Hill  Book  Company,  1969),  p.  19. 
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In  continuum  problems  of  any  dimensions,  the  field  variables 
(whether  pressure,  temperature,  displacement,  stress,  or  any  other  quantity) 
possess  an  infinite  number  of  values  because  it  is  a  function  of  each  generic 
point  in  the  body  or  solution  region.  Therefore,  the  problem  has  an  infinite 
number  of  unknowns.  The  finite  element  discredzation  procedure  reduces  the 
problem  to  one  of  a  finite  number  of  unknowns  by  dividing  the  solution  region 
into  elements  and  by  expressing  the  unknown  field  variable  in  terms  of 
assumed  approximate  functions  within  each  element. 

The  approximate  functions,  also  called  interpolation  functions  or 
shape  functions,  are  defined  in  terms  of  the  field  variable  values  at  specified 
points  called  nodes  that  usually  lie  on  element  boundaries  where  adjacent 
elements  are  considered  to  be  connected.8  In  addition  to  boundary  nodes,  an 
element  also  may  have  few  interior  nodes.  The  behavior  of  the  field  variables 
within  the  elements  is  completely  defined  by  the  nodal  values  of  the  field 
variables  and  shr.pe  functions.  For  the  finite  element  representation  of  a 
problem,  the  nodal  values  of  the  field  variable  become  the  new  unknowns. 
Once  these  unknowns  are  found,  the  shape  functions  define  the  field  variables 
through  the  assemblage  elements.9 


8K.  H,  Huebner,  The  Finite  Element  Method  for  Engineers  (London:  Wiley  and  Sons  Ltd.,  1975), 
p.  38. 

9Alghatam,  “An  Investigation  into  the  Use  of  the  Finite  Element  Method  for  Certain  Continuum 
Problems,”  pp.  25-32. 
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The  nature  of  the  solution  and  the  degree  of  the  approximation 
depend  not  only  on  the  size  and  number  of  the  elements  used,  but  also  on  the 
shape  functions  selected.  Functions  cannot  be  chosen  arbitrarily  because 
certain  compatibility  conditions  should  be  satisfied.  Often  functions  are  chosen 
so  that  the  field  variable  or  its  derivatives  are  continuous  across  adjoining 
element  boundaries. 

An  important  feature  of  the  finite  eloroent  method,  which  sets  it 
apart  from  other  approximate  methods,  is  the  ability  to  formulate  solutions  for 
individual  elements  before  putting  them  together  to  represent  the  entire 
problem.  This  means  that  complex  problems  can  be  reduced  by  considering  a 
series  of  greatly  simplified  problems. 

Another  advantage  of  the  finite  element  method  is  the  variety  of 
ways  in  which  one  can  formulate  the  properties  of  individual  elements  and 
even  change  them  within  the  solution  iterations.  This  applies  to  the  fluid  or 
continuum  properties  too. 

The  major  advantages  of  the  finite  element  method  over  the  better 
established  finite  difference  method  lies  in  the  ease  that  problems  involving 
complicated  irregular  geometries  and  arbitrary  local  mesh  refinement  can  be 
incorporated. 

The  present  finite  element  solution  of  the  three-dimensional 
Cartesian  coordinate  laser-induced  damage  to  mechanically  loaded  laminated 
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composite  panel  of  a  helicopter  is  compared  with  experimental  results  obtained 
from  the  NRL,  Washington,  D.C. 

It  was  desirable  in  this  thesis  to  have  the  finite  element  numerical 
method  of  solution  predict  the  physical  behavior  and  the  stress  patterns  and 
associated  heat  transfer  processes.  Great  emphasis  during  the  present 
investigation  has  been  placed  on  the  development  of  both  computational  and 
modeling  techmques  that  contribute  to  the  efficiency  of  the  finite  element 
method  as  applied  to  real-time  coupled  stress  and  heat  transfer.  A  major 
objective  of  this  work  is  to  show  the  applicability  and  flexibility  of  the  finite 
element  method  in  solving  continuum,  real-time  problems.  A  good  insight  into 
the  theory  and  applications  of  the  finite  element  method  can  be  found  in 
scientific  writings10,11’12,13  of  diverse  fields  of  engineering  such  as: 
vibration,  stress  analysis,  heat  transfer,  fluid  flow,  and  solar  energy.  In  order 
to  get  a  reasonable  grasp  of  the  foundation  knowledge  and  fundamentals  of  the 


10Mohammed  J.  Alghatam,  “The  Forced  Vibration  of  Beams  and  Frameworks  Using  the  Finite 
Element  Method”  (unpublished  BSC  Thesis,  Trent  Polytechnic,  Nottingham,  England,  UK,  1976),  pp. 
1-108. 

“Alghatam,  “An  Investigation  into  the  Use  of  the  Finite  Element  Method  for  Certain  Continuum 
Problems,"  pp.  1-178. 

12Mohammed  J.  Alghatam,  “Solar  Ventilation  and  Air-Conditioning  System  Investigation  Using 
the  Finite  Element  Method”  (Unpublished  PhD  Thesis,  University  of  Technology,  Loughborough,  Leics, 
England,  UK,  1986),  pp.  1-722. 

l3Mohammed  J.  K.  Alghatam,  “Solar-Powered  Air  Conditioning  System  Investigation  Using  Finite 
Element  Method,”  Solar  And  Wind  Technology.  Vol.  4,  No.  3  (Oxford:  Pergamon  Press,  1987),  pp.  243- 
268, 
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numerical  analysis  and  the  finite  element  method  and  techniques,  several 
specific  references14’16,16’17  provide  an  outstanding  source. 

The  four  basic  approaches  for  formulating  the  finite  element  to 
obtain  element  properties,18  are  as  follows: 

1.  The  direct  approach  (traceable  in  origin  to  stiffness  method  of 
structural  analysis); 

2.  The  variational  approach  (more  versatile  and  advanced); 

3.  The  weighted  residual  approach  (even  more  versatile  and 
advanced,  but  might  be  less  respectable  mathematically; 

4.  The  energy  balance  approach  (relies  on  the  balance  of  thermal 
and/or  mechanical  energy  of  a  system). 

Regardless  of  the  approach  used  to  find  the  element  properties,  the 
solution  of  a  continuum  problem  by  the  finite  element  method  is  always  in  the 
following  six  steps:19 


US.  Scheid,  “Numerical  Analysis,”  Schaum's  Outline  Series  (New  York:  McGraw-Hill,  1968)  pp. 
40-72. 

15J.  Timothy  Oden,  Finite  Elements  of  Non-linear  Continua  (New  York:  McGraw-Hill,  1972),  pp. 
1-131. 

18L.  J.  Segerlind,  Applied  Finite  Element  Analysis  (London:  Wiley  and  Sons  Ltd.,  1976),  pp.  12-98. 

17E.  Hinton  and  D.  R.  Owen,  Finite  Element  Programming  (Swansea:  Academic  Press,  1977),  pp. 
1-120. 

l8Alghatam,  “An  Investigation  into  the  Use  of  the  Finite  Element  Method  for  Certain  Continuum 
Problems,”  pp.  2-4. 

19Ibid.  pp.  2-4. 
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1.  Discretization  of  the  continuum:  dividing  the  domain  into 
elements; 

2.  Selection  of  interpolation  functions:  shape  functions; 

3.  Finding  the  element  properties:  matrix  equations  expressing 
element  properties  are  found  using  one  of  the  four  approaches 
just  mentioned; 

4.  Assembling  element  properties  to  obtain  the  system  equations: 
combining  matrix  equations  of  the  elements’  behavior  to  form 
overall  matrix  equations  expressing  the  behavior  of  the  entire 
solution  region,  domain,  or  system.  The  matrix  equations  for 
the  system  have  the  same  form  as  the  equations  for  an 
individual  element  except  that  they  contain  many  more  terms 
because  they  include  all  nodes.  Also  at  this  stage  the  boundary 
conditions  are  introduced  into  the  global  assembly  matrix.  The 
basis  of  the  assembly  is  that  at  shared  node  between  two  or 
more  elements,  value  of  the  field  variable  is  the  same  for  each 
element  sharing  the  node;20 

5.  Solution  of  the  system  equations:  solving  the  set  of  linear  or 
non-linear  simultaneous  partial  differential  equations; 

6.  Making  additional  computatioxis  if  desired:  sometimes  the 
solution  of  the  system  equations  is  used  to  calculate  other 

“Huebner,  “The  Finite  Element  Method  for  Engineers,”  pp.  49-56. 
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important  parameters,  e.g.,  calculating  heat  flux  from  known 
temperature  field. 

In  this  work  the  discretization  process  of  the  stress  part  is  done  via 
the  strain  energy  method  while  that  of  the  thermal  part  is  carried  out  by 
means  of  the  error  distribution  principle  (Galerkin  method),  using  weighted 
residuals,  where  the  weighting  functions  used  need  not  be  the  same  as  the 
interpolation  (shape)  functions  for  the  element  (Petrov-Galerkin).  In  each  case 
integration  by  parts  (Green’s  theorem)  is  used  to  invoke  the  boundary 
conditions.  The  governing  non-linear  partial  differential  equations  are 
transformed  into  “finite  element  equations”  governing  all  isolated  elements. 
These  local  elements  are  finally  collected  together  to  form  a  global  system  of 
non-linear  partial  differential  equations  or  “algebraic  equations”  with  proper 
boundary  and/or  initial  conditions  imposed.  The  use  of  suitable  boundary 
conditions  makes  the  coefficient  matrix  non-singular  and  not  ill-conditioned, 
thus  allowing  a  solution  to  be  obtained.  The  type  of  element  used  is  the  eight 
noded  isoparametric  quadrilateral  three-dimensional  (brick)  element 
(Serendipity  family).  Within  the  next  few  months  a  specially  formulated 
laminated  composite  finite  element  will  be  available  in  the  Finite 
Element/Personal  Computer  (F.E./P.C.)  package,  COSMOS/M  of  the  Structural 
Research  and  Analysis  Corporation  (SRAC),  used  in  this  thesis. 

This  thesis  consists  of  seven  chapters.  The  first  chapter  introduces 
the  subject  at  hand.  It  explains  the  scope  of  this  work  and  its  related  subjects. 


The  power  of  simulations  and  predictions  are  discussed  with  their  advantages 
and  limitations.  A  good  insight  into  the  finite  element  method  is  given  in 
conjunction  with  real  life  applications  of  engineering  continuum  problems. 

The  second  chapter  is  the  review  of  the  literature,  where  the  author 
conducted  a  thorough  research  of  available  unclassified  material  on  the  thesis 
subject.  It  introduces  the  few  leading  useful  military  sponsored  researches. 
It  discusses  their  accomplishments,  their  advantages,  limitations,  and  their 
relation  to  this  work. 

The  third  chapter  deals  with  the  structural  loading  boundary 
conditions  and  material  properties.  The  geometric  dimensions  of  the  laminated 
composite  panel  and  its  loading  boundary  conditions  are  configured  and 
explained.  The  experimental  critical,  buckling  and  failure  loads  from  NRL  are 
given  here.  The  composite  stress  analysis  in  laminated  composites  is  explained 
together  with  a  sketch  of  the  ply  orientation  and  staking  of  the  laminate 
matrix.  Also  given  is  the  ply  lay-up  and  thicknesses,  together  with  the 
material  mechanical  properties. 

The  fourth  chapter  explains  the  thermal  boundary  conditions  where 
the  laser  beam  geometry  and  energy  are  specified,  together  with  the  thermal 
properties.  The  material  degradation  thermal  relations  of  strength,  modulus, 
density,  specific  heat,  and  conductivity  are  shown  in  graphical  forms. 

The  fifth  chapter  describes  the  modeling  and  the  finite  element 
formulation  of  the  problem.  The  initial  considerations  and  the  development 


-  12  - 


plan  of  attack  for  the  solution  of  the  problem  is  given  there.  The  basic 
fundamentals  of  the  finite  element  method  as  applied  to  the  model  for 
subdivision  and  transformation  are  given.  The  three-dimensional  finite 
element  model  of  the  problem  is  drawn  and  shown.  The  theory  of  the 
composite  shell  element  is  introduced  together  with  its  elastic  constants, 
stress,  strain,  stiffness,  buckling,  and  load/strain  relationship.  Also,  in  this 
chapter  the  Galerkin  finite  element  formulation  of  the  general  energy  equation 
is  given  considering  conduction,  convection,  and  radiation.  The  solutions  were 
shown  for  linear  steady  state,  linear  transient,  and  non-linear  transient. 
Equilibrium  iterations  to  approach  the  correct  solution  were  explained  together 
with  the  convergence  criteria. 

The  sixth  chapter  presents  the  analysis  and  results.  Here  all  the 
relevant  input  and  output  data  is  displayed  for  both  the  stress  buckling  modes 
and  the  heat  transfer  analysis  together  with  the  subsequent  thermal  loads, 
contours,  and  other  parameters  are  shown.  They  are  indicative,  successful, 
accurate,  and,  perhaps,  impressive. 

The  seventh  chapter  contains  the  discussion,  conclusions,  and 
suggestions  for  further  work.  The  discussion  explains  the  analysis  and 
validates  the  results.  The  conclusions  are  summarized  in  accordance  with  the 
flow  of  the  chapters.  The  suggestions  for  further  work  indicate  the  possibilities 
of  improvements,  other  applications  to  be  tackled,  and  experimentation 
requirements. 
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CHAPTER  2 


REVIEW  OF  THE  LITERATURE 


The  author  has  conducted  deep  and  thorough  computerized  research 
of  the  existing  unclassified  literature  on  this  thesis  subject  in  the  United 
States,  Europe,  and  the  Soviet  Union.  However,  since  published  knowledge  in 
both  Europe  and  the  Soviet  Union  seems  to  be  primitive,  only  the  leading  few 
useful  military  sponsored  U.S.  research  documents  are  reviewed  here.  These 
documents  are  primarily  concerned  with  laser-induced  damage  experiments 
and/or  simulations  to  aircraft  or  missiles,  in  conjunction  with/without 
laminated  composites,  and  vice  versa.  This  chapter  discusses  those 
publications  accomplishments,  advantages,  shortfalls,  limitations  and  their 
relation  to  the  work  contained  in  this  thesis. 

Sheryl  K.  Bryan,  a  United  States  Air  Force  (USAF)  Aerospace 
Engineering  officer,  conducted  research  that  lead  to  her  degree  of  Master  of 
Science  in  1983.21  The  research  subject  was  a  reanalysis  method  for 

slSheryl  K.  Bryan,  “Reanalysis  Methods  for  Structures  with  Laser-Induced  Damage"  (Unpublished 
MS  Thesis,  Air  Force  Inst,  of  Tech.,  Wright-Patterson  AFB,  Ohio,  1983),  pp.  1-103. 


-  14  - 


structures  with  laser-induced  damage.  Strength  analysis  of  a  structure  which 
has  encountered  a  laser  strike  was  developed.  The  method  accounted  for  three 
types  of  laser-induced  damage: 

1.  Loss  of  structure  due  to  melting, 

2.  Change  of  material  properties  due  to  temperature  changes,  and 

3.  Addition  of  load  due  to  thermal  stress. 

The  program  used  heat  balance  calculations  over  successive  finite 
time  increments  on  an  array  of  finite  elements  bisecting  the  laser  beam  spot 
to  determine  the  temperature  distribution.  Those  results  were  then  converted 
to  structural  stiffness  parameters  and  the  structural  analysis  was  performed 
using  a  finite  element  based  reanalysis  method,  which  predicted  the  damage 
effects  from  the  initial  undamaged  solution.  The  researcher  claimed  that  the 
program  gave  agreement  with  results  obtained  from  a  separate  analysis.  No 
comparison  of  the  program  results  with  experimental  data  was  done.  Material 
used  was  isotropic  aircraft  metal  and  not  a  composite.  Buckling  was  not 
considered  and  the  heat  conduction  solution  was  two-dimensional. 

Laser- target  interaction  and  blast  wave  formation  in  DNA/NRL  laser 
experiment  was  numerically  simulated  by  a  team  of  researchers  in  1986.aa 
They  presented  details  of  a  one-dimensional  numerical  simulation  code 
intended  to  model  the  laser  equipment  in  high  pressure.  The  one-dimensional 

““U.S.  Department  of  the  Navy,  Numerical  Simulation  of  the  Laser-Target  Interaction  and  Blast 
Wave  Formation  in  the  DNA/NRL  Laser  Experiment,  by  John  L.  Giuliani  and  Margaret  Mulbrandon, 
Memorandum  Report  5762  (Washington,  D.C.:  Geophysical  and  Plasma  Dynamics  Branch,  Plasma 
Physics  Division,  Naval  Research  Laboratory,  1986),  pp.  1-56. 


-  15  - 


code  cannot  model  instabilities  development,  self-generated  magnetic  fields,  or 
hysteresis  losses.  The  results  indicated  unrealistically  large  amounts  of 
radiation  loss. 

A  high-order  theory  for  geometrically  non-linear  analysis  of 
composite  laminates  was  presented  in  1987.38  The  two  researchers  developed 
a  refined,  third  order  laminate  theory  that  accounted  for  transverse  shear 
strains.  Their  theory  allowed  parabolic  description  of  the  transverse  shear 
stresses,  and  hence  the  shear  correction  factors  of  the  usual  shear  deformation 
theoiy  were  not  required.  This  theory  also  accounted  for  small  strains  with 
moderately  large  displacements  (i.e.  Von  Karmen  strains).  They  derived 
closed-form  solutions  of  the  linear  theory  for  certain  cross-ply  and  angle-ply 
plates  and  cross-ply  shells.  The  finite  element  model  was  based  on 
independent  approximations  of  the  displacements  and  bending  moments  (i.e. 
mixed  formulations),  and  therefore  only  zero  continuity  approximations  were 
required.  The  mixed  variational  formulations  they  developed  suggested  that 
the  bending  moments  could  be  interpolated  using  discontinuous 
approximations  (across  inter-element  boundaries).  They  used  the  finite 
element  method  in  order  to  analyze  cross-ply  and  angle-ply  laminated  plates 
and  shells  for  non-thermal  stress  bending  and  for  natural  vibration.  This 


“J.  N.  Reddy  and  C.  F.  Liu,  A  Higher-Order  Theory  for  Geometrically  Nonlinear  Analysis  of 
Composite  Laminate.  NASA  Contract  Report  4056  (Blacksburg,  VA:  Virginia  Polytechnic  Institute  and 
State  University,  1987),  pp.  1-99. 
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theory,  however,  fell  short  of  catering  for  non-linear  material  or  thermal 
loading. 


Elastic  wave  dispersion  in  laminated  composite  plate  was  studied  by 
a  group  of  researchers  in  1987.w  They  examined  the  dynamic  behavior  of 
infinite  periodically  laminated  medium  using  the  stiffness  method  that  they 
have  developed  previously.  They  concentrated  on  studying  dispersion  of  waves 
in  a  laminated  plate.  In  their  approach  each  lamina  was  divided  into  several 
sublayers  and  the  displacement  distribution  through  the  thickness  of  each 
sublayer  was  approximated  by  polynomial  interpolation  functions  in  such  a 
way  that  displacements  and  tractions  were  continuous  across  the  interfaces 
between  adjacent  sublayers.  This  study  was  motivated  by  the  desire  to  model 
wave  propagation  in  a  continuous  fiber-reinforced  laminated  plate.  The 
deduction  was  that  if  the  wavelength  is  long  compared  to  the  fiber  diameters 
and  spacing,  then  each  lamina  can  be  modeled  as  a  homogeneous  transversely 
isotropic  medium  with  the  symmetry  axis  parallel  to  the  fibers.  The  overall 
effective  elastic  properties  of  such  a  medium  could  be  calculated  from  the  fiber 
and  matrix  properties  by  using  an  effective  modulus  theory  developed  in  1964. 
Such  an  assumption  has  been  made  in  their  paper.  Thus  the  dispersion 
characteristics  of  guided  waves  in  a  layered  anisotropic  plate  was  studied  by 
them.  In  the  analysis  they  assumed  that  each  lamina  is  transversely  isotropic 
with  the  symmetry  axis  aligned  with  either  the  x-axis  or  the  y-axis.  This 

a4S.  K.  Datta,  et  al.  Elastic  Wave  Dispersion  in  Laminated  Composite  Plate.  Contract 
N00014-86-K-0280,  (Bolder:  University  of  Colorado,  1987),  pp.  1*8. 


assumption  was  not  necessary  for  the  development  of  the  equations,  but  it  was 
made  to  keep  the  algebra  as  simple  as  possible  for  the  anisotropic  problem  they 
had.  The  results  gave  good  approximation  to  experimental  values.  Neither  a 
static  nor  a  thermal  loading  were  used. 

Funded  by  the  Department  of  Energy,  \  ashington,  D.C.,  Mr. 
Vincent  Prantil,  a  researcher  in  the  structural  mechanics  division  of  Sandia 
Labs.,  Livermore,  conducted  a  research  on  the  response  of  a  thin  cylindrical 
shell  under  lateral  impulse  loads  in  1988. 26  In  support  of  investigating  the 
lethality  levels  of  short  pulse  width  lasers,  this  combined 
analytical/experimental  program  was  underway  in  Sandia  National 
Laboratories  to  predict  large  deformation  response  of  conceptual  ballistic 
missile  designs  to  impulsive  loading.  Both  unpressurized  and  pressurized 
cylindrical  shell  models  have  been  analyzed  and  tested  The  type  of 
impulsively  loaded  metal  structures  in  question  exhibit  large  scale  plastic 
deformation  and  local  dynamic  pulse  buckling  instability  which  subsequently 
influence  the  loaded  surface  collapse.  Impulse  loads  were  simulated  using  both 
lead  spray  and  light  initiated  high  explosive  techniques.  Finite  element 
calculation  using  a  modified  form  of  the  Hughes-Liu  shell  element  based  on  a 
modified  code  of  the  DYNA3D  package  verified  several  characteristics  of  the 
structural  response  observed  by  the  researcher  in  experiments.  The  role  of  the 
computational  results  in  identifying  lethal  impulse  levels  and  response  nodes 

26V.  C.  Prantil,  Response  of  a  Thin  Cylindrical  Shell  Under  Lateral  Impulse  Loads  (Livermore: 

Structural  Mechanics  Division,  Sandia  National  Laboratories,  1988),  pp.  1-75. 
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were  described.  In  addition,  the  requirements  for  initial  imperfections  and  fine 
mesh  discretization  for  those  problems  were  discussed. 

The  laser  weapon/target  interaction  exhibited  three  characteristic 
response  time  regimes  which  were  treated  separately.  The  characterization  of 
the  full  response  was  broken  down  into: 

1.  Surface  physics  response, 

2.  The  one-dimensional  shock  wave  response, 

3.  The  overall  global  structural  response. 

The  surface  physics  addressed  the  relationship  between  weapon 
output  and  the  compressive  shock  wave  imparted  to  the  target  surface.  In  an 
x-ray  laser  encounter,  the  weapon  output  impinged  on  the  target  surface.  The 
radiant  energy  absorption  often  resulted  in  melting  of  the  exposed  target  skin 
and  plasma  blow  off,  delivering  a  compressive  shock  wave  into  target.  This 
exchange  took  place  in  a  submicrosecond  time  span.  The  response  of  the 
one-dimensional  shock  wave,  (also  called  material  or  hydrodynamic  response) 
was  concerned  with  several  effects.  Most  important  of  these  effects  was  the 
possibility  of  the  material  spallation  due  to  the  interaction  of  stress  waves  in 
the  material.  There  was  also  the  possibility  of  material  degradation  due  to  the 
short  intense  compressive  shock  wave  introduced  during  load  delivery.  The 
response  was  characterized  by  the  time  over  which  the  resulting  stress  waves 
propagate  through  the  shell  aluminum  alloy  metal  wall,  imparting  a  residual 
momentum  in  the  target.  Such  effects  occur  on  microsecond  time  scales  for  the 


target  materials  of  interest.  Finally,  the  structural  response  described  the 
overall  target  deformation  occurring  over  times  of  several  milliseconds  and 
beyond.  Within  this  last  time  frame,  the  entire  structure  is  responding  to  the 
applied  load.  This  is  exactly  what  we  are  trying  to  analyze  in  this  thesis,  but 
for  laminated  composites  instead  of  aluminum  alloys.  The  researcher 
investigates  and  describes  his  analytical  effort  to  capture  the  essential  aspects 
of  the  structural  response  modes.  In  doing  so  probable  relevant  failure 
mechanisms  could  have  been  identified.  Also  the  computational  capability  of 
non-linear  large  deformation  finite  element  codes  could  have  been  verified,  in 
addition  to  predicting  the  observed  behaviors  in  the  accompanying 
experiments.  Predicting  material  failure,  however,  has  not  been  attempted  by 
the  researcher,  and  actual  use  of  laser  was  not  experimentally  attempted. 

Failure  of  mechanically  loaded  laminated  composites  subjected  to 
intense  localized  heating  was  examined  by  Dr.  Ronald  Gularte  and  his 
research  team  in  the  NRL,  Washington,  D.C.,  in  1987.26  They  presented  a 
numerical  procedure  for  predicting  failure  in  laminated  composites  subjected 
to  simultaneous  intense  localized  heating  and  applied  mechanical  loads.  Their 
method  consisted  of  a  non-linear,  axisymmetric,  two-dimensional  finite 
difference  thermal  analysis  including  the  effects  of  surface  ablation,  as  well  as 
the  radiation  losses,  convection  losses,  and  temperature  dependent 


“Ronald  C.  Gularte  et  al..  “Failure  of  Mechanically  Loaded  Laminated  Composites  Subjected  to 
Intense  Localized  Heating,”  International  Journal  for  Numerical  Methods  in  Engineering.  Vol.  25 
(London:  Wiley  and  Sons  Ltd.,  1988),  pp.  561-570. 
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thermophysical  properties.  They  used  their  results  in  conjunction  with  a  flat 
plate  finite  element  code  to  predict  the  failure  of  laminated  composites. 
Predictions  based  on  the  numerical  analysis  were  compared  with  experimental 
results  obtained  from  mechanically  loaded  graphite  epoxy  coupons 
spot-irradiated  at  various  intensities. 

The  approach  presented  consisted  of  the  integration  of  a  high 
temperature  material  property  data  base  with  non-linear  thermal  analysis 
capability  to  produce  a  spatial  and  temporal  temperature  distribution  within 
an  intensely  locally  heated  structure  which,  in  turn,  was  used  in  the 
mechanical  code  to  predict  structural  failure.  The  thermal  degradation  of 
material  properties  resulted  in  material  non-linearities.  We  conclude  that  the 
overall  strategy  adopted  in  the  above  mentioned  research  to  be  considered  as 
an  uncoupled  iterative  procedure  between  a  stress  model  and  a  thermal  model. 
The  model  used  was  two-dimensional,  and  the  experimental  irradiated 
specimen  used  was  tensile  rather  than  compression  loaded.  There  were  some 
agreements,  but,  there  were  some  disparities  in  the  analysis  which  the 
researchers  attributed  to  the  slight  inaccuracies  in  the  thermal  calculations 
which  may  have  lead  to  errors  in  defining  the  temperature-dependent 
mechanical  properties.  They  added,  that  the  use  of  classical  plate  theory  to 
model  the  complex  stresses  state  in  the  vicinity  of  degraded  material  could 
have  been  an  over-simplification. 
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An  analytical  model  based  on  the  sublaminate  approach  and  fracture 
mechanics  was  developed27  to  study  the  growth  of  delamination.  This  report 
analyzed  the  inter-laminar  fracture  in  composites  under  combined  loading. 
Plane  strain  conditions  were  assumed  and  estimates  were  provided  for  the 
total  strain  energy  release  rate  and  modes  contributions.  The  energy  release 
rate  estimated  were  used  in  combination  with  experimental  data  on 
graphite/epoxy  laminates.  Reasonable  agreement  was  demonstrate 1  for  the 
range  were  the  experimental  observations  indicated  transverse  crack  tip 
delamination  to  be  the  predominant  failure  mode.  The  analytical  model  was 
two-dimensional,  and  the  assumptions  of  plane  strain  conditions  with 
negligible  thickness  strain  are  over-simplification.  No  thermal  analysis  was 
done. 

A  final  technical  report  describing  the  research  done  on  damage 
models  of  continuous  fiber  composites  was  introduced  by  Mr.  D.H.  Allen, 
Aerospace  Engineering  Department,  Texas  A&M  University  in  1988.28  This 
report  was  submitted  to  the  Air  Force  Office  of  Scientific  Research,  USAF.  The 
report  summarized  research  completed  during  a  four  year  period  under  USAF 
grant.  The  objective  of  the  research  has  been  to  develop  an  accurate  damage 
model  for  predicting  strength  and  stiffness  of  continuous  fiber  laminated 


^Erian  A  Armanios,  Analysis  of  Interlaminar  Fracture  in  Composites  Under  Combined  Loading 
(Atlanta:  Georgia  Institute  of  Technology,  1988),  pp.  1-30. 

“David  H.  Allen,  Research  on  Damage  Models  for  Continuous  Fiber  Composites.  Final  Technical 
Report  (College  Station,  Texas  A&M,  1988),  pp.  1-86. 
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composite  media  subjected  to  fatigue  or  monotonic  loading  and  to  verify  this 
model  with  experimental  results  obtained  from  composite  specimens  of  selected 
geometry  and  make-up.  The  work  performed  was  in  four  parts: 

1.  Development  of  the  constitutive  equations  relating  stresses  to 
strains  and  damage  internal  state  variables  which  might  have 
been  used  in  a  stress  gradient  field, 

2.  Development  of  internal  states  variables  growth  laws  as  a 
function  of  load  history  for  matrix  cracking,  interlaminar 
fracture,  etc., 

3.  Development  of  finite  element  algorithms  capable  of  evaluating 
ply  properties  in  damaged  components, 

4.  Perform  experiments  on  components  with  selected  stacking 
sequences  in  order  to  verify  the  model. 

On  the  basis  of  the  experimental  studies,  a  general  framework  was 
developed  using  continuum  damage  mechanics  to  characterize  the  response  of 
laminates  with  matrix  cracks.  This  model  was  compared  favorably  with 
experimental  evidence.  The  model  was  then  extended  to  account  for  both 
matrix  cracks  and  delaminations,  and  compared  favorably  with  experimental 
results.  The  analyses  were  all  two-dimensional,  and  no  thermal  effects  or 
loading  were  introduced.  This  research  concentrated  on  the  micro-structural 
damage  rather  than  the  overall  or  macro-  structural  damage  which  is  of 
engineering  design  interest  to  us. 
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An  interim  report  on  the  development  of  a  progressive  failure  model 
for  strength  of  laminated  composite  structure  was  composed  by  P.Y.  Tang,  for 
the  Naval  Ocean  Systems  Center,  San  Diego  in  1989.39  The  researcher 
deduced  that  the  various  progressive  failure  schemes/models  results  do  not 
compare  satisfactorily  with  experimental  data  when  geometrically  complicated 
laminated  composite  specimens  are  in  question.  The  researcher  proceeded 
with  an  assumption  based  on  the  introduction  of  an  improved  ply  failure 
criterion  into  the  progressive  failure  model  should  improve  the  accuracy  of 
predicting  the  strength  of  laminated  composite  structure.  The  improved 
version  of  the  finite  element  package  COSMIC/NASTRAN  of  NASA  is  proposed 
to  be  used,  in  conjunction  with  a  Stanford  University  package  called  PDHOLE 
which  was  specially  written  to  predict  the  tensile  strength  of  a  composite 
laminate  strip  containing  an  open  hole.  This  proposal  by  the  above  mentioned 
researcher  is  to  take  place  in  the  fiscal  years  1990-91  timeframe  to  try  and 
prove  the  researcher  assumption.  No  thermal  loading  were  considered,  and  the 
proposed  modeling  is  in  two-dimensions. 

Compression  failure  mechanics  in  unidirectional  composites  was 
studied  in  a  NASA  technical  memorandum.30  Possible  failure  modes  of 


29P.  Y.  Tang,  Development  of  a  Progressive  Failure  Model  for  Strength  of  Laminated  Composite 
Structure  (San  Diego:  Naval  Oceans  System  Center,  1989),  p.  2. 

30H.  Thomas  Hahn  and  Jerry  G.  Williams,  Compression  Failure  Mechanics  in  Unidirectional 
Composites.  NASA  Technical  Memorandum  85834  (Langley,  VA:  NASA  1984),  pp.  1-26. 
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constituent  material  were  summarized  and  analytical  models  were  reviewed. 
The  composites  used  were  unidirectional  and  no  heating  was  used. 

The  buckling  modes  and  failure  characteristics  of  differently 
stiffened  composite  shear  panels  were  described  in  NASA  Technical 
Memorandum  2153,31  Shear  stress-strain  at  tow  temperatures  were 
compared  with  linear  buckling  analysis  and  failure  behavior  is  discussed. 


31Mark  J.  Stuart  and  Jane  A.  Hagman,  Buckling  and  Failure  Characteristics  of  Graphite-Polvimide 
*  Shear  Panels.  NASA  Technical  Memorandum  2153  (Langley,  VA:  NASA,  1983),  pp.  1-21. 
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CHAPTER  3 


STRUCTURAL  LOADING  BOUNDARY 
CONDITIONS  AND  MATERIAL  PROPERTIES 


In  this  chapter  the  structural  loading  boundary  conditions  are 
inscribed  and  material  properties  are  stated.  The  geometric  dimensions  of  the 
laminated  composite  panel  and  its  loading  boundary  conditions  are  configured 
and  explained.  The  experimental  critical,  buckling  and  failure  loads  from  the 
NRL  are  given  here.  The  composite  stress  analysis  in  laminated  composites 
is  explained  together  with  a  sketch  of  the  ply  orientation  and  staking  of  the 
laminate  matrix.  Also  given  is  the  ply  lay-up  and  thicknesses,  together  with 
the  material  mechanical  properties. 

As  we  discussed  in  the  introduction  we  would  like  to  predict  the 
buckling  and  post-buckling  behavior  of  heated  flat  laminated  composite 
rectangular  orthotropic  compression  panels  with  edges  rotationally  restricted. 
A  good  example  of  one  particular  application  for  such  a  high  temperature, 
low-weight  material  is  the  aft  body  of  the  Space  Shuttle. 


-26- 


COMPOSITE  STRESS  ANALYSIS: 


A  lamina  (ply),  as  shown  in  figure  3.2(a),  is  a  flat,  sometimes  curved 


as  in  a  shell,  arrangement  of  unidirectional  fibers  or  woven  fibers  in  a  matrix. 


A  laminate,  as  shown  in  figure  3.2(b),  is  a  stack  of  laminae  with  various 
orientations  of  fiber  directions  called  ply  orientations  in  the  laminae.32 


PLY 

ORIENTATION 

ANGLE 

(degree) 


Figure  3.2  Laminated  Composite 


Due  to  the  inherent  heterogeneous  nature  of  fiber-reinforced 
composite  materials,  they  are  studied  in  both  micromechanics  and 
macromechanics.  In  the  first,  composite  material  behavior  is  studied  by 
examining  the  interaction  of  the  constituent  fibers  and  matrixes  on  a 
microscopic  scale.  In  macromechanics,  composite  material  behavior  is  studied 

3aP.  Y.  Tang,  Development  of  a  Progressive  Failure  Model  for  Strength  of  Laminated  Composite 
Structure,  p.  2. 
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by  presuming  the  material  to  be  homogeneous  and  detecting  the  effect  of  the 
constituent  material  to  meet  a  particular  structural  requirement.33  This 
latter  study  is  what  is  important  to  us  in  this  work,  from  the  engineering 
design  point  of  view. 


“Ibid,  p.  2. 
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MATERIAL  PROPERTIES: 


Graphite/Epoxy 
Hercules  AS4/3501-6 


MECHANICAL  PROPERTIES 

Units 

Tvoe  AS4 
Fiber 

AS4/3601-6 

Comnosite 

Tensile  Strength 

Ksi 

520 

315 

0  Degree  Direction 

MPa 

3,587 

2,170 

Tensile  Modulus 

Msi 

34 

21 

0  Degree  Direction 

GPa 

235 

145 

Short  Beam  Shear: 

Strenrtti 

Ksi 

_ 

18.5 

Unidirectional 

Mpa 

128 

Flexural  Strength 

Ksi 

260 

Unidirectional 

Mpa 

1,793 

Compression  Strength  at  room  temperature  is:  64  Ksi 


INTER-LAMINAR  SHEAR: 

At  Room  Temperature  7.3  Ksi 

At  200  Degrees  F  8.8  Ksi 
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GEOMETRIC  CONSIDERATIONS: 


In  order  to  reduce  the  computing  time,  effort  and  cost,  only  a  quarter 
of  the  laminated  composite  panel  was  considered  to  be  modeled  (upper  right 
section)  due  to  symmetry.  In  addition,  the  material  properties  were  such  that 
the  quadrant  or  jc/2  modeling  is  more  than  adequate  within  the  engineering 
bounds  of  the  study,  as  opposed  to  detailed  complete  circle/panel  modeling. 

The  uniformity  and  circumference  accuracy  of  the  laser  beam 
incident  circle  on  the  panel  is  adequate  for  engineering  design  purposes. 
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CHAPTER  4 


THERMAL  BOUNDARY  CONDITIONS 

The  thermal  boundary  conditions  together  with  the  laser  beam 
geometry  and  energy  are  specified  in  this  chapter.  Also  stated  are  the  material 
thermal  properties.  The  material  degradation  thermal  relations  of  strength, 
modulus,  density,  specific  heat  and  conductivity  are  presented  graphically. 

CASE  SPOT  SIZE/RADIUS  (in) 

1  1.000 

2  2.000 

(*)3  3.000 
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Laser  Beam  Radius 


Figure  4.1  Irradiation  Boundary  Conditions 


* 


... ^  The  degradation  of  the  room  temperature  properties  with  increasing 

temperature  is  given  in  the  following  next  five  figures/graphs.  The  first  graph 
shows  Ultimate  Strength  (%  R.T.  value)  against  Temperature  (Deg.  R),  while 
the  second  graph  shows  Modulus  (%  R.T.  values)  against  Temperature  (Deg. 
R),  the  third,  fourth,  and  fifth  graphs  show  the  temperature-dependent  density, 
specific  heat,  and  conductivity  of  graphite/epoxy.34 


^Ronald  C.  Gularte,  Head,  Engineering  Materials  Group,  Naval  Research  Lab.,  and  G.  Holderby, 
Wright  Research  and  Development  Center,  Washington,  D.C.,  private  communications,  1990-1991. 
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THERMAL  PROPERTIES  UNITS  TYPE  AS4 


AS4/3501-6 


Fiber  Composite 


Specific  Heat,  C 
(RT) 

Thermal  Conductivity 
0  Pea,  direction 
(RT) 

90  Deg,  direction 
(RT) 


J/G.Deg.K 
Btu/h.ft.Deg.F 
or:  Cal /g  Deg.C 

Btu/h.ft.Deg.F 
Watt/M  Deg.K 

Btu/h.ft.Deg.F 
Watt/M  Deg.K 


0.92 

3.1 

0.22 

0.74 

15 

2.65 

26 

4.60 

0.36 

mmm 

0.62 

We  are  going  to  limit  ourselves  to  one  Loading  and  One  Heating 
Condition  as  indicated  by  (*)  in  text  in  Chapters  (3,4). 

Nx  (Critical)  =  2410  lb/in,  Nx  (Buckling)  =  3927  lb/in 
A  (6M)  Spot  Diameter  with  an  Irradiance  of  (1.000  KW/Sq.cm) 
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THERMAL  CONSIDERATIONS: 


Both  ablation  and  spall  at  the  interfaces  of  the  laser  beam  circle  on 
the  panel  were  neglected  due  to  the  fact  that  the  temperature  generated  in  this 
exercise  does  not  bring  the  panel  to  these  stages.  This  was  set  by  the  (NRL), 
however,  it  is  possible  to  include  the  effect  of  both  readily  in  the  program 
should  this  be  required.  The  instantaneous  damage  that  takes  place  in  the 
panel  is  due  to  mechanical  buckling  rath  ^r  than  anything  else.  If  elevated 
temperature  causing  spalling  and  ablation  are  reached  then  the  effect  of  both 
would  be  secondary,  since  the  structural  buckling  damage  would  have  already 
took  place.35 

The  effect  of  air  flow  and  painting  on  the  surface  to  resemble  real 
life  helicopter  outer  panel  situation  were  both  examined.  The  effect  of  both  are 
trivial  due  to  the  speed  of  the  heating  process.36,37 


“Bryan,  “Reanalysis  Methods  for  Structures  with  Laser-Induced  Damage,”  pp.  48-70. 
“Ibid,  pp.  48-70. 

37Prantil,  “Response  of  a  Thin  Cylindrical  Shell  Under  Lateral  Impulse  Loads,”  pp.  51-58. 
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PROPERTIES  DEGRADATION 
WITH  INCREASED  TEMPERATURE 


8TRk8S  TYPE 
SHEAR 

~+~  LONQTUONAL  C0MPRE88N 
TRAN8VER8S  TENSION 
-B-  TRANSVERSE  COMPRE8SN 
LONQTUONAL  TENSION 


ULTIMATE  •TRINVTH  (%  R.T.  VU.UI) 


TEMPERATURE  (Deg.  R) 
LASER  INDUCED  DAMAGE 


Figure  4.2  Properties  Degradation  (1) 


PROPERTIES  DEGRADATION 
WITH  INCREASED  TEMPERATURE 


0  600  1000  1500  2000  2600  3000 


TEMPERATURE  (Deg.  R) 

8TRE38  TYPE 

TRAN8VER8E  AND  SHEAR  —  LONGITUDINAL 

LASER  INDUCED  DAMAGE 

Figure  4.3  Properties  Degradation  (2) 
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TEMPERATURE  DEPENDENT  DENSITY 
OF  GRAPHITE/EPOXY  LAMINATED  COMPOSITE 


g/cublc  cm 

:  1 

. 

- 

Density 

i  i  i  i  i 

O  500  1000  1500  2000  2600  3000  3500 

TEMPERATURE  (Deg.  C) 

THERMAL  RESPONSE  TO  LASER  HEATING 

Figure  4.4  Temperature  Dependent  Density  of  Composites 


TEMPERATURE  DEPENDENT  SPECIFIC  HEAT 
OF  GRAPHITE/EPOXY  LAMINATED  COMPOSITES 


J/g.Deg  C 

8 1 - 


Specific. heat 

Epoxy  DaandMtlQD  CaDacitv  (Co) 

-  s 

(343-5101  DaaC 

_ - —  Or  aphis  Subllmat loo 

Ta  •  3316  Daa  C 

.  - 

Ha f  43  KJ/g 

1 - 1 _ L . 1 _ 1 - 1 - 

0  600  1000  1600  2000  2600  3000  3600 

TEMPERATURE  (Deg.  C) 

THERMAL  RESPONSE  TO  LASER  HEATING 

Figure  4.5  Temperature  Dependent  Specific  Heat  of  Composites 
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TEMP.  DEPENDENT  THERMAL  CONDUCTIVITY 
OF  GRAPHITE/EPOXY  LAMINATED  COMPOSITE 


0  600  1000  1600  2000  2600  3000  3600 

TEMPERATURE  (Deg.  C) 

—  Trentveree  Direction  — Fiber  Direction 

THERMAL  RESPONSE  TO  LASER  HEATING 

Figure  4.6  Temp.  Dependent  Thermal  Conductivity  of  Composites 
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CHAPTERS 


MODELING  AND  THE  FINITE  ELEMENT 
FORMULATION  OF  THIS  PROBLEM 

The  modeling  and  the  finite  element  formulation  of  this  thesis 
problem  are  presented  in  this  chapter.  The  initial  considerations  and  the 
development  plan  of  attack  for  the  solution  of  the  problem  are  given  here.  The 
basic  fundamentals  of  the  finite  element  method  as  applied  to  the  model 
subdivision  and  transformation  is  given.  The  three-dimensional  finite  element 
model  of  the  problem  is  drawn  and  shown.  The  theory  of  the  composite  shell 
element  is  introduced  together  with  its  elastic  constants,  stress,  strain, 
stiffness,  buckling,  and  load/strain  relationships.  Also,  in  this  chapter  the 
Galerkin  finite  element  formulation  of  the  general  energy  equation  is  given 
considering  conduction,  convection  and  radiation.  The  solutions  were  shown 
for  linear  steady  state,  linear  transient,  and  non-linear  transient.  Equilibrium 
iterations  to  approach  the  correct  solution  were  explained  together  with  the 
convergence  criteria. 

In  order  to  be  able  to  capture  all  the  influencing  behavior  in  this 
complicated  set-up  of  both  stress  and  temperature  distribution  dynamic 
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dependence  on  material  properties,  the  author  chose  to  work  in 
three-dimensions  for  both  the  stress  and  thermal  analysis.  The 
three-dimensional  finite  element  mesh  used  to  model  the  current  problem  is 
shown  on  the  next  page.  Since  the  geometrical  configuration  of  the  laminated 
composite  panel  in  question  is  symmetrical  in  all  three  cartesian  coordinate 
axis,  obviously,  only  quarter  of  the  panel  is  modeled.  However,  the  initially 
suggested38  development  plan  of  attack  for  the  solution  of  the  problem  was 
only  in  two-dimensions  and  is  as  follows: 

1.  To  conduct  an  analysis  to  verify  that  Nx-Critical  is  indeed  the 
critical  buckling  load. 

2.  Perform  a  simplified  (axi-symmetrical  or  one-dimemional) 
thermal  analysis  which  should  result  in  the  temperature 
distribution  as  a  function  of  time. 

3.  Finally,  use  the  thermal  results  as  input  to  the  indpiently 
critically  loaded  model  and  determine  the  amount  of  time 
required  to  initiate  buckling. 

4.  If,  time  allows,  one  could  use  the  previous  approach  on  the 
post-buckled  panel  to  predict  the  time  required  to  fail  the 
irradiated  panel. 

The  basic  fundamental  of  the  finite  element  method  is  the 
sub-division  of  the  domain  of  interest  into  smaller  bounded  subdomains  called 


“Gularte  and  Holderby,  Private  communications. 
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elements.  Such  elements  are  connected  together  at  discrete  points  on  their 
extremities  called  nodes. 

The  behavior  of  a  dependent  variable  over  each  element  may  be 
approximated  by  continuous  functions  such  that  at  points  X,Y,Z: 

■  £t"W.  (5,1a) 

where  n  is  the  number  of  elemental  nodes.  Similarly  derivatives  may  be 
expressed,  typically: 


f£(WZ)  -  £  ^ 

C/A  OA 


dNt 


(5.1b) 


In  the  finite  element  method  N(X,Y,Z)  are  referred  to  as  the  shape 
(or  interpolation)  functions.  These  are  generally  cast  in  a  parametric  form 
using  an  (£,r|,0  coordinate  system,39  and  in  continuum  mechanics  are  chosen 
to  ensure  functional  continuity  at  element  boundaries.  However,  such 
continuity  need  not  exist  for  function  derivatives,  these  elements  are  referred 
to  as  (Co)  continuous. 

The  transformation  of  derivative  terms  from  parametric  to  the 
physical  coordinate  system  is  obtained  by  the  chain  rule  and  by  invoking 
equation  (5.1b),  typically: 


dX  A 

as  "  h  as  ' 


38Zienkiewicz,  “The  Finite  Element  Method,”  pp.  132-256. 


-  42  - 


Then  adopting  matrix  notation,  in  two  dimensions: 


dNt 

ax 

& 

dNt 

35 

a? 

ax 

dNt 

ax 

ar 

dNt 

.dr\. 

.3n 

fri. 

.ar 

i.e. 


then 


dNt 

dNt 

-m 

ax 

dNt 

dNt 

.an. 

IF 

dNt' 

3W, 

ax 

-  \J\~l 

35 

dNt 

dNt 

[ar, 

3ti 

(5.2) 


When  the  dependent  variable  is  interpolated  at  the  same  nodes  as 
the  dependent  variables  (X,Y,Z  coordinates)  the  element  is  said  to  be 
isoparametric.  When  the  dependent  variable  is  interpolated  at  fewer  nodes, 
the  element  is  described  as  superparametric.  Some  typical  element  geometries 
are  shown  in  Figure  6.1.  Further  information  concerning  the  elements  type, 
shape  functions,  continuity,  compatibility,  etc.,  can  be  found  in  “The  Finite 
Element  Method.”40 


40Ibid.  pp.  132-256. 
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The  three-dimensional  finite  element  model  of  the  problem  is  shown 
in  the  next  Figure  5.2.  The  convention  for  stacking  sequence  in  local  element 
coordinate  and  the  temperature  distribution  across  the  thickness  of  the  plate 
is  shown  for  the  composite  shell  element  in  Figure  5.3.  The  side  view  of  the 
model  used  is  comprised  of  120  shell  elements  for  each  layer  of  the  laminated 
composite  panel  quadrant. 


Figure  5.1  Element  Types 

(a)  2D  Linear,  (b)  2D  Quadratic,  (c)  3D  Linear,  (d)  3D  Quadratic 
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Figure  5.2  Three-dimensional  Finite  Element  Model 


5.1  COMPOSITE  SHELL  ELEMENT: 


Figure  5.3  Composite  Shell  Element 
Convention  for  stacking  sequence  in  local  element  coordinate  and 
the  temperature  distribution  across  the  thickness  of  the  plate 


5.: 


ELASTIC  CONSTANTS: 


Although  for  the  most  general  cases  of  anisotropy  the  number 
of  independent  elastic  constants  is  21,  these  constants  are  required  to  define 
the  general  anisotropic  material  in  pounds  per  square  inch.  This  number  is 
considerably  reduced  if  the  material  internal  composition  possesses  any  kind 
of  symmetry.  A  state  of  anisotropy  possesses  three  mutually  orthagonal  planes 
of  symmetry  at  each  layer.  If  the  reference  system  of  orthagonal  axes  (X,Y,Z) 
is  parallel  to  the  principle  material  axes  (1,2,3),  and  (z  =  3)  axis  is 
perpendicular  to  the  midplane  of  the  structure  at  each  point,  the  elasticity 
matrix  D  relating  stress  and  strains  is  obtained  as: 

Dl  Du  0  0  Olfe, 

d12  d2  0  0  0  62 

0  0  D,  0  0  y12 

0  0  0  D4  0  Y 13 

0  0  0  0  D5  Y23 


(5.3) 


where  y  equals  shear  stress,  G  equals  shear  modulus,  and 
Di  -  £i/d-v12v21) 

^2  ”  ^j/(l  “^12 ^21^ 

&12  "  -V12V2l) 

D%  “  G12 
Da  -  Kx  x  Gu 

Ds  -  XjxGjj 
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the  terms  Kj  and  K,  are  shear  correction  factors  which  are  chosen  to  match  the 
plate  theory  with  certain  classical  solutions.41 

If  the  principle  axes  of  anisotropy  1,  2  do  not  coincide  with  the 
reference  axes  X,  Y,  but  are  rotated  by  a  certain  angle  theta  (0),  as  illustrated 
in  Figure  5.4  below: 


Z  Y 


x  -  rj  +  rj  +  r3k 
y  -  s,»  +  s2j  +  s3k 
Z  -  bxi  +  bj  +  b3k 

C  -  Cosd  -  rxbx  +  r2b2  +  r3b3 
S  -  Sind  -  b1sl  +  b2s2  +  b3s3 


I,  J ,  k  are  unit  vectors  used 
to  define  a  vector  quantity  in 
toms  of  its  components. 

hlf  b2,  b3  are  components  of 
the  z  vector  as  slf  Sj,  S3  and 
ri»  r2>  h  are  components  of 
the  y  and  x  vectors  respectively. 


Figure  5.4  Principle  Axes  of  Anisotropy 
Sign  convention  for  principle  in-plane  material 
axes  (1,2)  and  local  element  axes  (X,Y) 


4lS.  P.  Timoshenko  and  J.  M.  Gere,  Theory  of  Elastic  Stability.  2nd  Ed.  (New  York:  McGraw-Hill, 
1961),  pp.  418-419. 
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The  new  elasticity  matrix  D  is  determined  by  the  following 


transformation. 


°1.2.3 

ei,2,3 


To 


*,y.i 


r'e. 


*.7.1 


(5.4) 


where 


ax,y,t  "  {  °r«  V  *«>  }* 

e*.w  -  { e*’  v  y*.  y*  }r 
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7/  o 
0  t2 


and 


c2 

S2 

SC 

52 

C2 

-sc 

-2SC 

2SC  C2-S2 

c  s 

o 

_ 1 

t2  - 

r'  - 

-s  c 

[0  t3 j 

o  -  r_l  o 

*.7.* 

1.2.3 

"  rl2?eU.3 

-  T^DJV. 

WL 

0  ■ 

*./.*  n  . 

-  Z)e 

*.7.* 

0 

(5.5) 


(5.6) 
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where  [QJ  is  the  bending  part  of  the  material  matrix  and  [QJ 
relates  to  the  shear  effects. 


5.1.2  FORCE  RESULTANT-STRAIN  RELATIONSHIP: 


The  stress  in  the  K*  layer  in  terms  of  the  reference  surface 
strains  and  curvature  are: 


• 

r 

'  ex«  ' 

°y 

V 

•  +  z ' 

JC, 

► 

.  T*r . 

V. 

V 

(5.7) 


T 


n 


Stress  resultants  and  stress  couples  can  now  be  defined  as  those  for 
homogeneous  plate  by  integrating  the  above  equations  over  each  layer  and 
summing  the  resulting  expressions  over  the  n  plate  layers. 


zk*i 


N,’Ny>Nxf  "  E  /  (°*>W)<Jt)  dz 


k-l  z. 


(5.8) 
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-  r  /  e>,.<VV*  zdz 


*- 1  Z 


Zfl 


v*r,  ’  £  /  <*«.**>”  * 

*-i  zt 
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where  Z  equals  the  distance  from  the  reference  surface;  Nx>  Ny>  and  N*y 
represent  the  in-plane  shear  forces;  M*,  My,  depict  the  in-plane  normal  and 
shear  moments;  and  Vx,  Vy,  and  portray  the  transverse  shear  forces. 
Carrying  out  the  indicated  integrations  yields  equation  5.9  listed  below: 


iln  A13  1  Bn  Bu  B13  I 

I  1 

N, 

^ 21  -^22  ^*23  i  &21  Bn  0 

il31  j432  j433  |  B31  Bn  B33  | 

Mx 

i  j 

^11  fi13  1  Al  ^12  -®13  1 

M, 

-®21  -®22  ^23  j  D2l  ^22  ^23  j  0 

M 

■*3' 

"®31  “®32  ^33  j  ^31  ^32  ^33  | 

1  1 

|  |  £11  *12 

»i  J 

o  1  0  1  *n  *»  . 

where 


(5.9) 


^  -  E  «?*„>'  Cz*.i  -  z.) 

*-l 

B*  -  ^E«?vr  <z»*>  -  z«> 

Ak-\ 


(5.10) 


The  presence  of  [B]  above  implies  coupling  between  bending  and 
extension.  This  term  is  equal  to  zero  if  the  laminate  is  symmetric  about  the 
reference  surface. 


5.1.3  LOCAL  ELEMENT  STIFFNESS  MATRIX: 


The  strain  energy  of  the  plate  is  given  by: 


« -  iff 


e0' 

r 

ABO 

•  •> 

e 

it 

- 

B  D  0 

i 

k 

.  Y  . 

0  0  E 

Y  , 

dxdy 


(5.11) 


which  can  be  written  as: 


’■iff 


,0  )  T 


A  B 
B  D 


-?U}T 


*  |  tedy  +  ±ff  {  Y  }T  {  £  H  Y  }  dxdy 


A  B 


//[*»  f  [B>]dxdy  +  [l[B,]T[E][Bt]dxdy 


(5.12) 

\{U) 


where 


(5.13) 


-51- 


m  i  o  ■ 
---I — 

[*.]- 
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o  |  [as] 

.  0  1  [flfll 

6x18 

1 

(5.14) 


and  [BP],  [BB],  and  [BS]  are  the  kinematic  matrices. 

[BP]  and  [BB]  are  partitioned  matrices  used  for  bending  and  [BS]  is 
the  same  for  shear. 


The  stiffness  matrix  is  given  by: 


[K\  -  [Kb]  +  [&] 


(5.15) 


where 

m  ffm7  [Bb}<h<fy 

i 

[&]  -  //[ib]r[£] 


5.2  BUCKLING  AND  TN-PLANE  FORCES: 

When  a  structure  is  subjected  to  in-plane  forces  (i.e.  forces  that 
act  tangent  to  a  member  axis  or  mid  surface),  membrane  strain  energy 
develops  inside  the  structure.  The  effect  of  the  membrane  strain  energy  on  the 


-  52  - 


total  stiffness  of  structure  is  considered  by  adding  the  geometric  stiffness 
matrix  to  the  global  stiffhess  matrix. 


The  geometric  stiffhess  matrix  is  independent  of  elastic 
properties  and  depends  only  on  element’s  geometry,  displacement  held  and 
state  of  stress.  It  accounts  for  the  effect  of  existing  forces  on  the  bending 
stiffhess  of  structure. 

In  presence  of  membrane  forces,  the  equation  of  static 

equilibrium: 


M  M  - « 


will  change  to: 

[K  +  Ka]  {£/}  -  {*} 


(5.16) 


(5.17) 


where  [K]  is  the  stiffhess  matrix,  [Kg]  is  the  geometric  stiffhess 
matrix  for  the  given  load  vector  {R},  and  {U}  is  the  displacement  vector. 

As  can  be  seen  from  the  graph  in  Figure  5.6,  representing  the  beam 
of  Figure  5.5,  the  transverse  displacement  becomes  very  large  and  approaches 
infinity  as  the  axial  load  approaches  a  certain  limit.  This  limit,  known  as 
Pcriud>  is  the  buckling  load  for  the  structure. 
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Figure  5.5  •  Cantilever  Beam 

Buckling  occurs  when  a  member  of  the  structure  converts  membrane 
strain  energy  into  strain  energy  bending.  After  buckling,  the  member  will 
usually  fail  because  sudden  large  bending  deformations  are  needed  to  store 
membrane  energy.  If  the  same  problem  is  solved  in  the  absence  of  the 
transverse  loading  shown  in  Figure  5.6,  the  load-displacement  curve  is 
obtained  as  in  Figure  5.7. 

P 


Load-displacement  Curve  for  a  Typical  Beam/Column 
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Axial  Force 

P 


Figure  5.7 

Load-displacement  Curve  for  a  “Perfect”  Column 


Note  that  in  this  case  of  Figure  5.7  above,  the  is  a  bifurcation 
point,  after  which  two  possible  equilibrium  configurations  exist.  In  the 
post-buckling  state,  path  AC  of  the  above  Figure  5.7  is  never  followed,  because 
it  is  the  path  of  greatest  resistance. 

In  order  to  formulate  the  buckling  problem,  equation  5.17  can  be 
rewritten  in  the  form  of: 


(K  +  Ajy  {<?}  -  {R} 


(5.18) 


{Q}  is  the  transformed  {U}  vector  needed  to  convert  equation  5.17 
into  equation  5.18. 
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where  is  the  geometric  stiffness  matrix  for  an  arbitrary  state  of  stress  in 
the  structure  and  X  is  an  orbitary  scaler  multiplier. 


If  the  above  system  is  in  equilibrium,  then  addition  of  a  virtual 
displacement  pattern,  (U)  would  keep  the  two  configurations  in  equilibrium, 
(i.e.): 

( [jq  ♦  [K0]  )  [U]  -  {R} 

(5.19a) 

or 


(K  +  X[jq  )  {U}  -  0 


(5.19b) 


Equation  5.19b  has  the  form  of  a  general  eigenvalue  problem  where 
X  is  equal  to  the  load  multiplier  necessary  to  create  buckling  and  {U}  is 
buckling  mode  shape.  It  should  be  noted  that  the  vector  {U}  is  only  a  mode 
shape  which  gives  the  general  pattern  of  displacement,  and  that  the  value  of 
{UJ  at  a  certain  point  has  no  physical  significance  by  itself.  Also  it  is 
important  to  notice  that  the  solution  of  equation  5.19b  gives  the  overall 
buckling  load  of  the  structure  as  per  this  thesis  requirements,  and  buckling  of 
individual  members  cannot  be  predicted  by  this  ^  ation.  A  negative  value  for 
X  signifies  that  the  overall  structure  is  in  tension  and  that  the  direction  of  all 
loads  acting  on  the  structure  must  be  reversed  to  cause  buckling. 


5.3  gale: 


RKIN  FINITE  ELEMENT  FORMULATION  OF  THE 
HEAT  CONDUCTION  (GENERAL  ENERGY)  EQUATION; 


The  governing  partial  differential  equation  of  heat  conduction 
(general  energy)  in  the  cartesian  coordinates  is  given  by: 


_a_ 

dx 


+  Qg  -  pc 


36 

dt 


(5.20) 


The  associated  boundary  conditions  are: 

1.  Dirichlet  boundary  conditions  (prescribed  temperatures)  0  =  8* 

on  boundary  SI.  (5.21) 

2.  Neuman  boundary  conditions  (prescribed  heat  flow) 


*>•30  «■  30  p-  30  0  n 

-  *•  n"* '  'V'  ‘  Q ' 


3 z 


(5.22) 


on  boundary  S2. 


The  Dirichlet  boundary  conditions  in  the  present  analysis  are  treated 
as  a  special  case  of  Neuman  boundary  conditions  by  assuming  a  large  value  of 
convective  heat  transfer  coefficient  and  an  ambient  temperature  which  is  equal 
to  the  prescribed  temperature. 


Note  that  the  value  Q  in  the  Neuman  boundary  conditions 
represents  the  total  heat  flux  conducting  into  the  body  through  the  boundary 
on  which  Neuman  boundary  conditions  are  specified.  For  a  simple  case  with 
no  phase  change  condition,  the  conservation  of  heat  flow  at  the  boundary  can 
be  stated  as: 

Total  heat  flux  input  -  Total  heat  flux  conducting  into  the  body 
The  total  heat  flux  conducting  into  the  body  may  be  expressed  as: 

Q-Qe  +  Qr  +  Qp 

where 

Qe  -  Total  heat  flux  due  to  convection 

-  K(0e  -  0) 

Qr  -  Total  heat  flux  due  to  radiation 

-  A,(0„  -  6) 

where  the  equivalent  radiation  coefficient  is  given  by: 

h,  -  Be  (e*,  -  e2)(0,  +  B) 

Qp  -  Total  prescribed  heat  flux 


(5.23) 

(5.24) 

(5.25) 

(5.26) 
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Divide  the  domain  into  N  finite  elements.  In  each  of  these 
finite  elements,  assume  the  variation  of  temperature  6  as  a  function  of 
interpolations  Hj  and  the  nodal  temperatures: 


-  [  H(x,y,z)  ]8W 


(5.27) 


In  Galerkin  method,  the  integral  of  weighted  residue  over  the 
domain  is  set  to  zero.  The  weights  in  this  method  are  the  same  as  the 
interpolation  functions  where  ///  means  a  volume  integral,  and 
//  means  fn  a  surface  integral. 


dxdydz 


0 


(5.28) 


Applying  Green-Gauss  (integrating  by  parts)  theorem  and  using  the  equations 
(5.22  to  5.26),  the  equation  (5.28)  can  be  rewritten  in  the  following  form: 


d  . .  «,  30  9  , .  „  901  .  .  . 


(5.29) 


.  Iff  HfiC  dxdydz  +  / [  H,  (A,  *  kr)6  ds 


‘  / f!  Hfi,  +  //  h<  (M.  *  M,) 


The  above  system  of  equations  (for  i  =  1,  N)  can  be  expressed  in 
matrix  form  as: 


/ / [«]  ^ + *,)[«]<&]  a  *  [fjfpcmm**  ] 


39 

at 


-  [///mrc,*<H  *  [//mr(Acfle**,e,)&]  ♦  [  //[fl]r 9,* 


(5.30) 


Thus  the  finite  element  system  of  equations  can  be  rearranged  in  the 
following  form: 


[  [*»]  +  PJ  *  [*r]  ]  9  *  W  *  -  9,  *  9,  *  9.  ♦  ?, 


[*j  3  +  [C]  S  -  q 


(5.31) 


Heat  conduction  matrix 


[*»]  '///  [fi]r  [Kt]  [3]  dxdydz 

(5.32) 


Thermal  conductivity  matrix 


Thermal  gradient  matrix 


Thermal  Load  Vectors 


Internal  heat  generation  load  vector 


«•„  -  f f f  [H]T  Q,  dxdydz 


(5.38) 


Convective  heat  transfer  load  vector 


f'  -  //  K  W  * 


(5.39) 


Radiation  heat  transfer  load  vector 

f,  -  //  K  W  a,  * 


(5.40) 


Prescribed  heat  flow  load  vector 

(5.41) 

5.3.2  SOLUTION  PROCEDURES 

5.3.2.1  LINEAR  STEADY  STATE  ANALYSIS: 


Since  the  temperature  is  not  a  function  of  time 

|  (8)  -  ° 

*  (5.42) 


Equation  (5.31)  can  be  expressed  as 


[  [**]  ♦  fc]  ]  ®  +  <?c  +  1P 

(5.43) 
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Note  that  the  radiation  boundary  conditions  are  not  present  in  the 
above  equation  because  this  is  a  linear  analysis.  Also  note  that  the  material 
properties  and  boundary  conditions  are  not  functions  of  either  time  or 
temperature. 

Steady  state  nodal  temperature  distribution  can  be  obtained  by 
solving  the  above  system  of  equations  after  applying  the  Dirichlet  boundary 
conditions.  The  Dirichlet  boundary  conditions  can  be  easily  applied  without 
disturbing  the  order  of  the  matrices  as  follows: 

Let  the  prescribed  temperature  of  node  “i”  be  0*.  The  two  steps 
procedure  is: 

1.  Replace  the  diagonal  element  of  ith  row  of  effective  conduction 

matrix  [K]  by  a  large  number: 


i.e.,  Ku  -  1.0U+15 


(5.44) 


2.  Replace  the  corresponding  thermal  load  by  the  product  of  this 
large  number  and  the  prescribed  temperature  i.e. 


qt  -  (l.OE+15)  6* 


(5.45) 
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Then  solve  the  finite  element  system  of  equations  to  get  the  steady  state 
temperature  distribution. 

S.3.2,2  LINEAR  TRANSIENT  ANALYSIS: 

Analysis  with  constant  material  properties  and  constant 
boundary  conditions  with  no  radiation  come  under  this  category.  Recalling  the 
finite  element  system  of  equations: 

[A]  0  +  [C\  0  -  4  (5.46) 

Here  0  is  the  time  derivative  of  the  temperature.  Assuming 
constant  gradient  in  the  neighborhood  of  time  t,  9  may  be  expressed  in  any 
of  the  following  methods: 

Forward  Difference  Method 


0  -  [  ('^>0  _  *0  ]/a* 
Backward  Difference  Method 


(5.47) 


0  -  [  ‘0  .  (»-a*)0  ]/a* 

Central  Difference  Method 


(5.48) 
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0.[('*A»)g_  ("AOg]/Af 


(5.49) 


In  general,  0  can  be  expressed  as  a  function  of  a  parameter  a  as 


follows: 


-  [  (»-*»>  0  -  «9]  /  (aAf) 


(5.50) 


since  material  properties  are  constant 


<'♦•*'>[*]  -  *[K\ 


<r^«A0[C]  -  f[q 


(5.51) 


(5.52) 


Now  consider  the  heat  flow  equilibrium  at  time,  t  +  ocAt, 


f+aA,)[K]  (,+aAO0  +  (f<^A»)0  .  (t+aAt)g  ^ 


Substituting  5.50,  5.51,  and  5.52  into  5.53  to  get: 


m +  [c]  e  -  (f+aA<)*  -  ^  [c]  * 


(5.54) 
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[£*]  (,+aA,)g  -  (**•**)£  +  tg 

(5.55) 

With  the  known  load  vector  at  time,  t  +  aAt,  the  temperature 
distribution  at  time,  t  +  aAt,  can  be  obtained  by  solving  the  equation  5.55. 

Then  the  temperature  distribution  at  time,  t  +  At,  can  be  obtained 
as: 

('+A0g  -  (1  -  iVs  +  1  C'+«A‘>g 

l  a)  a  (5.56) 

Note  that  matrix  [K*]  is  independent  of  the  temperature  distribution,  i.e.  [K*] 
is  constant  throughout  the  solution.  Hence,  the  triangularization  of  the  matrix 
[K*]  is  to  be  done  only  once  for  the  entire  solution  process. 

At  each  time  step,  the  load  vector  corresponding  to  the  time,  t  +  aAt, 
and  the  temperature  distribution  at  same  time  step,  t  +  aAt,  can  be  obtained 
by  back  substitution.  Then  the  temperature  at  time,  t  +  At,  can  be  obtained 
using  eq  ration  5.56.  A  step  by  step  procedu  re  to  solve  linear  transient  heat 
conduction  problems  is: 

Initial  calculations 


1.  Form  heat  conduction  matrix  [KJ 
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2.  Form  heat  convection  matrix  [KJ  (if  required) 


3.  Form  heat  capacity  matrix  [C] 


4.  Compute  the  effective  conduction  matrix  as: 

M  -  N +  W +  ^ 

5.  Modify  [K*]  for  temperature  boundary  conditions 


6.  Compute  the  initial  resistance  vector 


7.  Triangularize  effective  conduction  matrix  [K*] 

For  each  solution  time  step,  follow  the  solution  below: 

1.  Form  the  heat  flow  vector  at  time,  t  +  aAt 

2.  Compute  the  resistance  vector 

tf.-L-tq* 

o.  Compute  the  effective  heat  flow  vector 

+  tg 
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4.  Modify  load  vector  for  temperature  boundary  conditions 

5.  Solve  the  heat  flow  equilibrium  equations  for  the  temporature 
distribution  at  time,  t  +  ocAt: 

[K*]  (*+*A*>g  - 


6.  Compute  the  temperature  distribution  at  time,  t  +  At,  as: 

C'+A*>0  -  +  1 

7.  Return  back  and  proceed  to  step  1  to  find  the  temperature 
distribution  for  the  next  time  step.  Otherwise,  stop  the  solution 
procedure  according  to  convergence  criteria  or  procedure  ending 
concept. 


S.3.2.3  NON-LINEAR  TRANSIENT  ANALYSIS: 


equations: 


Recall  the  governing  system  of  the  finite  element 


[  [**]  *  [*.]  *  pq  ]S  *  [C]8  +  (5.57) 

Rearrange  the  above  system  of  equations  in  the  following  form: 
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Rearrange  the  above  system  of  equations  in  the  following  form: 


pqB  ♦  [q8  -  f 


where  the  equivalent  conduction  matrix  is: 

w  -  [«*]  ♦  [*.]  *  ft] 


and  the  equivalent  thermal  load  vector  is: 


4-  +  $C  +  +  tip 


(5.58) 


(5.59) 


(5.60) 


It  is  obvious  that  the  temperature  dependent  material  properties  and 
the  radiation  boundary  conditions  make  this  system  of  equations  highly 
non-linear.  Hence,  it  is  appropriate  to  use  incremental  solution  procedure,  i.e. 
calculation  of  temperature  distribution  at  time,  t+ocAt,  from  the  known 
temperature  distribution  at  time,  t.  Hence,  consider  the  heat  flow  equilibrium 
at  time,  t+ocAt: 


('+«A0[Jf]  (t+«Af)g  _  ( t*aAt)g 

where  the  left  superscript  denotes  the  time. 


(5.61) 
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Since  the  matrices  [K]  and  [C]  are  functions  of  temperature,  the 


solution  process  can  be  started  with  the  following  linearization: 

-  *[K\ 

<,+«A'>[q  -  *[C] 

('♦■AO0  „  <g  +  a§(P) 


(5.62) 

(5.63) 

(5.64) 


where  A0^  is  a  vector  of  incremental  temperatures  between 
times  t  and  t+ocAt.  The  linearized  system  of  equations  can  be  represented  as: 


(,+aA,>5  +  ‘[C]  (f+aA,)0  - 


(5.65) 


Assuming  linear  variation  of  temperature  over  the  time  increment: 


(*♦«*»>  0  -  'g 
aA  t 


(5.66) 


Substitute  5.66  into  5.65  and  rearrange  terms  to  get: 


w  *  ^  w 


5^(0)  _  (t+aAt)g  _  tg 


(5.67) 
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solution  of  this  system  of  equations  gives  the  first  approximation  of 
incremental  temperatures.  Hence,  the  approximate  temperatures  at  time, 
t  +  aAt,  are  given  by: 


(6.68) 


Therefore,  the  temperatures  at  time,  t  +  At,  are  computed  as: 

(t* A*)0<P)  _  1  +  ^  ±) 

«  {  (5.69) 

5.3.2.4  EQUILIBRIUM  ITERATIONS  TO  APPROACH  THE 
CORRECT  SOLUTION  (CONVERGENCE): 

Solution  of  the  equation  6.67  does  not  give  the  correct 
incremental  temperatures  because  of  the  linearization  of  equation  6.61,  i.e.  the 
substitution  of  this  solution  into  equation  5.61  gives  an  unbalanced  heat  flow. 
Hence,  this  solution  has  to  be  updated  to  take  care  of  this  unbalanced  heat 
flow.  This  procedure  should  be  repeated  until  the  residual  heat  flow  is  within 
the  required  tolerance.  These  iterations  which  achieve  the  heat  flow 
equilibrium  are  called  equilibrium  iterations. 

The  temperature  correction  at  the  end  of  iteration  “i”  may  be 


expressed  as: 


(f«Ai)g(0  .  (r+«A<)g(<-»)  + 


(5.70) 


and 


A le®  -  IS®  -  IS<,_1) 


(5.71) 


where  I 2e®  is  tlxe  correction  in  incremental  temperatures  during  iteration 
“i.”  Now,  consider  the  heat  flow  equilibrium  of  equation  5.61  at  the  end  of 
iteration  “i”: 


(r+«Af)[jq  (t+aA0g(0  + 


[  <»♦« Ar)g<0  _  rg 


a  At 


-  (f«Aj)^ 


(5.72) 


Substitute  equation  5.70  into  5.72  to  get: 


(r+«A0[jq  + 


1 

aAf 


„  (f+«A»)^  _  (f+«Ar)^j  (»+«Af)g(<_1)  _  (f+aAf)|CJ  (f+aArJg1 


(<-l) 


(5.73) 


Corrections  for  incremental  temperatures  may  be  obtained  by  solving 
the  equation  5.73.  On  close  observation,  it  is  obvious  that  this  method  is 
similar  to  the  Newton-Raphson  solution  of  equation  5.61.  Note  that  this 
method  requires  the  triangularization  of  the  coefficients  matrix  in  each 
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iteration  which  is  very  expensive.  Hence,  a  modified  Newton-Raphson  scheme 
is  followed  in  which  the  coefficients  matrix  is  not  updated  in  each  iteration. 
Thus  each  iteration  only  involves  a  back  substitution  which  is  very 
inexpensive. 

The  governing  equations  for  the  modified  Newton-Raphson  iterations 
are: 


w 


1 

«A  t 


'[CJ 


_  (f+«A»)^  _  (f+cAOrgj  (»*«A0g 


(<-l) 


(5.74) 


A#0  -  W'X)  +  AA^ 


(5.75) 


(r*«Af)g(0  _  (f+aA»)g(<-»)  +  440(0 


(5.76) 


Repeat  the  above  calculations  until  convergence  is  achieved. 

5.3.2.5  CONVERGENCE  CRITERIA: 


It  is  assumed  that  the  solution  has  converged  if  the 
following  criteria  is  satisfied: 
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(5.77) 


I  aa*o  L 

- -  <  Convergence  tolerance 

|  (t..A*>g(0|2 

where  the  Convergence  tolerence  is  0.00001  for  this  thesis,  and 
the  norm  of  vector  V  is: 

I  via  -  /FfTvfn fTTTi? 

(5.78) 

Vv  V2,  ...Vu  being  the  components  of  the  Vector  V.  The  temperature 
distribution  at  time,  t  a  At,  may  be  obtained  from  the  temperature  distribution 
at  time,  t  +  ccAt,  from  the  equation: 

< _  1  (t+aAt)fi(Q  +  .  JL'j 

«  l  “J  (5.79) 

Repeat  this  procedure  to  find  the  temperature  distribution  at  all  other  required 
time  steps.  A  step  by  step  procedure  to  solve  the  non-linear  transient  heat 
conduction  problems  is  as  follows: 

Initial  calculations 


1.  Define  initial  temperature  vector 


2.  Initialize  time  step  counter 
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3. 


Form  conduction  matrix  [KJ 


4.  Form  convection  matrix  [KJ  (if  required) 

5.  Form  radiation  matrix  [KJ  (if  required) 

6.  Form  heat  capacity  matrix  [C] 

7.  Compute  the  effective  conduction  matrix: 

[*1  -  FI  -  FI  *  FI  * 

8.  Modify  [K*]  for  temperature  boundary  conditions 

9.  Compute  the  initial  heat  flow  vector  due  to  conduction 

*.-['[*»]♦  'FI*  'Fj3 *s 

10.  Triangularize  the  matrix  [K*] 

For  each  solution  time  step,  follow  the  solution  procedure  given 
below: 


1.  Increment  the  time  step  counter 
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KSTEP  =  KSTEP  +  1 


2.  Increment  the  time 
TIME  =  TIME  +  aAt 

where  At  is  the  time  step  increment.  In  case  of  variable  time 
step  analysis,  substitute  the  appropriate  time  step  increment 
in  the  above  equation. 

3.  Check  for  matrix  reformation? 

NO:  Advance  to  step  9 

IF  (KSTEP.EQ.1)  advance  to  step  9 
YES:  Continue 

4.  Form  heat  conduction,  convection,  radiation  and  heat  capacity 
matrices  corresponding  to  the  temperature  distribution  at  the 
end  of  the  last  time  step. 

5.  Compute  the  effective  conduction  matrix  as: 

m  -  M +  w +  m  *  ir,[C] 

6.  Modify  fK*]  for  temperature  boundary  conditions 
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7.  Compute  the  internal  flow  vector  due  to  conduction 

*  -  { 'N  -  ♦  W  « 

8.  Triangularize  matrix  [K*] 

9.  Define  the  heat  flow  vector  at  time,  t  +  aAt,  i.e.  (  ) 

10.  Compute  the  effective  heat  flow  vector 

_  tp 

11.  Solve  the  following  sys.cm  of  equations  for  nodal  temperature 
increments  a3w  by  back  substitution 

12.  Check  for  heat  flow  equilibrium? 

YES:  Advance  to  step  16 

NO:  Continue 

13.  Find  the  temperature  distribution  at  time,  t  +  aAt 

(f+«A»)g  „  'g  +  Af)^ 


14.  Find  the  temperature  distribution  at  time,  t  +  At 


15.  Check  for  the  total  number  of  solution  time  steps. 

IF  (KSTEP.LT.NSTEP)  Go  to  step  1 
Otherwise,  stop  the  solution 

16.  Initialize  iteration  counter  1=0 

17.  Compute  the  temperature  distribution  at  the  beginning  of  the 
first  equilibrium  iteration 

(r+«A«)g(0)  _  tg  +  A0<o> 

18.  Increment  the  iteration  counter 
1  =  1+1 

19.  Compute  the  nodal  temperatures  and  the  time  derivatives  at 
the  end  of  the  previous  iteration 


20.  Calculate  the  out  of  balance  heat  flow  vector 

21.  Compute  the  effective  heat  flow  vector 

<f+*At)q*  .  (t+aLt)g  _ 

22.  Solve  for  the  corrections  in  temperature  increments  by  back 
substitution 

[jt]aa0«  -  (<+bA^* 

23.  Compute  the  total  increments  in  temperature  and  temperature 
at  time,  t  +  ocAt 

A8(0  -  aB(<_1)  +  AA0W 

m  (f+«Af)g(<-l)  +  ^0(0 

24.  Check  for  convergence 

|  AA0©  L 

- - -  <  Convergence  tolerance 

|  (i+«A»)0(O  | 


If  convergence:  advance  to  step  25 

IF  (I  .LT.  Max.  no.  of  iterations)  go  to  step  18 

Otherwise,  restart  using  smaller  time  step  or  increase  the 

maximum  allowable  number  of  iterations 

25.  Compute  the  temperature  distribution  at  the  end  of  time,  t  + 
At 

('♦*00  -  1  (f«Af)g(0  +  ^  _  _Lj  $ 

26.  Check  for  the  total  number  of  solution  time  steps? 

IF  (KSTEP.LT.NSTEP)  go  to  step  1 
Otherwise,  stop  the  solution. 
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CHAPTER  6 


ANALYSIS  AND  RESULTS 


This  is  the  chapter  of  analysis  and  results.  Here,  all  the  relevant 
input  and  output  data  are  displayed  for  both  the  stress  buckling  loads  and  the 
heat  transfer  analysis  together  with  the  subsequent  thermal  stress  analysis. 
Color  presentations  of  stress  and  thermal  modes,  contours,  and  other 
parameters  are  shown.  They  are  indicative,  successful,  accurate  and  perhaps 
impressive. 

Enclosed  are  the  relevant  input  and  output  data  from  the  input  and 
output  files  of  some  of  the  problems  run  by  the  COSMOS/M  program.  A  total 
of  thirteen  drawings  are  displayed.  The  program  uses  shell  of  revolution 
elements  to  solve  the  buckling  problems  subjected  to  symmetric  and 
asymmetric  static  and  dynamic  loads.  Significant  amounts  of  engineering  time 
are  saved  by  using  shell  of  revolution  elements  rather  than  general  shell 
elements.  The  external  load  is  decomposed  into  Fourier  components  and  the 


stress  solution  is  composed  of  these  Fourier  components  or  a  summation  of 
them.  Buckling  is  solved  in  terms  of  Fourier  components. 

The  geometric  modeler  part  defines  the  key  points  and  boundaries 
to  create  the  problem  3-D  mesh  interactively.  The  advance  dynamic  analysis 
performs  a  step  by  step  time  history  analysis,  frequency  response,  mode 
shapes,  and  buckling  response.  Heat  transfer  analysis  is  fully  done  for  the 
transient  state. 

BUCKLING: 

The  Figures  6.1,  6.2,  and  6.3  show  buckling  loads  and  modes  of  the 
laminated  composite  panel  as  had  been  described  in  the  original  drawings.  As 
mentioned  before,  cnly  a  quarter  of  the  panel  was  considered  (upper  right 
section)  due  to  symmetry.  Two  types  of  boundary  conditions,  namely  simply 
supported  and  fixed,  were  accounted  for.  The  buckling  loads  indicated  for  each 
case  are  in  terms  of  pounds  per  unit  length  (lbs/in).  For  the  (simply-supported) 
boundry  conditions  all  four  edges  were  considered  as  such,  with  three  edges  of 
the  panel  being  free  to  move  in  the  x-direction.  All  four  edges  were  considered 
fixed  for  the  (fixed-fixed)  boundry  conditions  except  for  motion  in  the  x- 
direction  along  three  edges. 
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THERMAL; 


The  heat  transfer  analysis  was  performed  by  applying  the  specified 
heat  flux  on  the  circular  region  in  the  middle  of  the  panel  in  accordance  with 
the  time  curve  specified  by  the  NRL.  The  thermal  output  results  were 
displayed  graphically  for  several  different  time  steps.  Figures  6.4  through  6.7 
show  the  temperature  distribution  for  time  steps  1, 10, 12,  and  25  respectively. 

Subsequently  a  thermal  stress  analysis  was  performed  by  applying 
the  temperature  distribution  acquired  from  the  heat  transfer  analysis.  Figures 
6.8  through  6.13  show  the  stress  distribution  for  temperatures  of  time  steps  1 
and  25.  Three  components  of  stress  are  considered  for  each  case. 

INPUT  AND  OUTPUT  FILES  (DATA  AND  RESULTS): 

1.  BUCK1SES:  An  input  file  giving  in  the  data  for  the  simply 
supported  boundary  conditions  of  the  laminated  composite  panel  in  static 
buckling. 

2.  BUCK2SES:  An  input  file  giving  in  the  data  for  the  encastr6 
(fixed-fixed)  boundary  conditions  of  the  laminated  composite  panel  in  static 
buckling. 
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Figure  6.1 

Simply-Supported  Boundary 
Condition  for  Buckling  Mode 


node 
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Figure  6.2 

Simply-Supported  Boundary 
Condition  for  Buckling  Mode 
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Figure  6.3 

Fixed-Fixed  Boundary  Condition 
for  Buckling  Mode  1 
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Figure  6.4 

Temperature  Di s tri buti on 
Time  Step  1 
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Figure  6.5 

Temperature  Distribution 
Time  Step  10 


nee 


89 


Figure  6.6 

Temperature  Distribution 
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Figure  6.7 

Temperature  Distribution 
Time  Step  25 


I 
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Figure  6.8 

Stress  in  X  Direction 
Time  Step  1 


Figure  6.9 

Stress  in  Y  Direction 
Time  Step  1 
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Figure  6.10 

Von  Mises  -  Time  Step  1 
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Figure  6.11 
Stress  in  X  Direction 
Time  Step  25 
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Figure  6.12 
Stress  in  Y  Direction 
Time  Step  25 
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Figure  6.13 


3.  Tl.LQQ:  An  input  file  giving  in  the  material  properties  and 
specified  temperature  and  flux  boundary  conditions  to  go  from  heat  transfer 
to  stress.  Such  properties  included  Young’s  Modulus  of  Elasticity  in  each 
cartesian  coordinate  direction,  Poisson’s  ratio  in  each  cartesian  coordinate 
direction,  thermal  conductivity  in  each  cartesian  coordinate  direction,  specific 
heat  capacity  in  each  cartesian  coordinate  direction,  material  density  in  each 
cartesian  coordinate  direction,  and  coefficient  of  heat  transfer  (Thermal 
Expansion)  in  each  cartesian  coordinate  direction. 

4.  BUCKl.OUT:  An  output  file  for  the  buckling  of  the  simply 
supported  boundary  conditions  of  the  laminated  composite  panel.  It  gives  out 
the  system  data  with  regard  to  the  number  of  equations  solved,  matrix 
elements,  half  bandwidth,  Max.  &  Min.  diagonal  stiflhess  matrix  values, 
number  of  eigenvalu  .  and  mode  shapes.  It  prints  out  for  each  mode  shape 
the  X,  Y,  and  Z  displacements  as  well  as  rotations.  At  the  end  of  this  output 
file  it  lists  the  time  consumed  in  performing  each  solution  operation. 

o.  BUCK2.QUT:  An  output  file  for  the  buckling  of  the  encastr6 
(fixed-fixed)  boundary  conditions  of  the  laminated  composite  panel.  It  gives 
out  the  system  data  with  regard  to  the  number  of  equations  solved,  matrix 
elements,  half  bandwidth,  Max.  &  Min.  diagonal  stiflhess  matrix  values, 
number  of  eigenvalues,  and  mode  shapes.  It  prints  out  for  each  mode  shape 
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the  X,  Y,  and  Z  displacements  as  well  a  rotations.  At  the  end  of  this  output 
file  it  lists  the  time  consumed  in  performing  each  solution  operation. 

6.  Tl.SES:  An  input  file  for  the  heat  transfer  session.  It  gives  in 
geometry  and  material  properties. 

7.  Tl.MOD:  An  input  file  for  the  heat  transfer  module.  It  gives 
in  the  material  properties,  prescribed  nodal  temperatures  from  the  previous 
heat  transfer  analysis  together  with  the  degrees  of  freedoms  and  the 
constraints. 


8.  Tl.TEM:  An  output  file  giving  the  temperature  results  from 
the  heat  transfer  analysis.  It  lists  the  modeling  mesh  characteristics,  solutions 
time,  and  it  lists  the  nodal  temperatures  for  each  time  step. 

9-  Tl.OUT:  An  output  file  giving  the  stress  output  at  each 
temperature/time  step,  i.e.  thermal  stresses  output.  This  file  lists  both  the 
principle  stress  (SIGMA)  and  the  shear  stress  (TAU)  in  each  of  the  cartesian 
coordinate  directions  (X,  Y,  and  Z)  for  each  node.  It  also  lists  out  the  axial  end 
rotational  force  moments  in  each  of  the  cartesian  coordinate  directions  (X,  Y, 
&  Z).  At  the  end  a  small  solution  times  log  is  displayed. 


-98- 


CHAPTER  7 


DISCUSSION,  CONCLUSIONS,  AND 
SUGGESTIONS  FOR  FURTHER  WORK 


This  seventh  chapter  is  the  discussion,  conclusions,  and  the 
suggestions  for  further  work.  The  discussion  follows  the  sequence  of  the 
chapters  and  discusses  the  analysis  explaining  and  validating  the  results.  The 
conclusions  are  summarized  in  accordance  with  the  flow  of  the  chapters.  The 
suggestions  for  further  work  indicate  the  possibilities  of  improvements,  other 
applications  to  be  tackled,  and  experimentation  requirements. 

DISCUSSION 

Laser  directed  energy  weapon  systems  are  relatively  a  new  field. 
The  design,  effects,  and  experience,  as  well  as  knowledge  of  operation,  are  not 
widespread.  However,  such  knowledge  and  expertise  can  be  found  in 
governmental  and  commercial  forms  although  the  literature  is  not 
comprehensive. 
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It  was  stated  at  the  beginning  of  this  thesis  that  one  of  the  main 
objectives  was  to  demonstrate  the  applicability  and  flexibility  of  the  finite 
element  method  in  solving  continuum  problems  in  applied  military 
engineering.  The  present  results  were  compared  with  those  obtained  from 
experimental  NHL  measurements.  The  good  agreement  proves  tho 
applicability  of  the  finite  element  technique  in  the  present  class  of  problems. 

The  finite  element  technique  incorporating  strain  energy  method  was 
used  to  solve  the  static  buckling  problem.  On  the  other  hand,  incorporating 
the  Galerkin  weighted  residual  approach,  and  primitive  variables  formulation, 
is  used  to  solve  the  variable  conductivity  general  energy  equation. 

Results  of  the  stress  analysis  showed  accurately  the  loading 
distribution,  while  the  temperature  results  showed  clearly  and  accurately  the 
temperature  gradients  affecting  the  thermal  stresses  leading  to  buckling  and 
consequent  failure.  Modeling  the  puise  buckling  response  is  critical  as  the 
structure  is  directly  and  instantaneously  weakened  leading  to  substantial 
collapse  of  the  mechanically  and  thermally  loaded  panel. 

From  Figures  6.1,  6.2,  and  6.3  it  can  be  seen  that  the  buckling  load 
for  the  simply  supported  (S.S.)  case  showed  a  half-wave  deformation  at  mode 
8  whereas  for  mode  1  which  gives  the  lowest  buckling  load  two  waves  appeared 
for  its  mode  shape.  Thus,  we  consider  the  buckling  mode  1  with  lowest  critical 
load  as  9.26  pounds  per  square  inch  for  the  S.S.  case  and  16.9  pounds  per 
square  inch  as  the  critical  load  for  the  fixed-fixed  case  in  the  first  mode  of  buckling. 
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The  simply  supported  case  emulates  a  wing  or  a  tail  of  a  helicopter 
while  the  fixed-fixed  case  models  a  part  of  the  main  body.  Since  each  time 
step  in  the  analysis  is  one  (1)  second  corresponding  to  temperature  steps  then 
buckling  critical  stress  is  reached  in  time  step  10  i.e.  ten  seconds  with  our 
1KW  per  square  centimeter  laser  load  making  the  panel  lowest  temperature 
of  625  F,  while  at  the  irradiation  section  of  the  panel  the  maximum 
temperature  is  3750  F. 

The  stress  analysis  was  performed  for  temperature  distributions  that 
were  obtained  from  the  transient  heat  transfer  analysis.  This  analysis  was 
made  fcr  temperatures  at  time  steps  1  and  25.  The  program  allows  the 
reading  of  temperature  quantities  at  different  time  steps  and  use  them  for 
thermal  stress  analysis.  The  stress  units  are  in  pounds  per  square  inch  and 
the  Von-Mises  stress  corresponds  to  an  average  equivalent  stress  used 
principally  for  design  purposes  and  calculated  in  terms  of  other  components. 

The  negative  stress  means  compressive  stress  and  the  positive  (+Ve) 
stress  means  tensile  stress.  The  thermal  stress  analysis  was  performed  for  the 
S.S.  model  for  temperature  conditions  at  time  steps  1  and  25,  since  we  are 
considering  firing  at  helicopter  tails.  However,  the  required  exposure  time  of 
10  seconds  is  relatively  long  for  weapon  system  active  tracking  while  firing  at 
a  moving  target  in  the  same  time.  The  simplest  solution  is  to  increase  the 
power  of  the  laser  beam  from  1KW  per  square  centimeter  to  5KW  per  square 
centimeter  which  will  result  in  a  panel  buckling  time  of  no  more  than  one 
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second  which  correspond  to  the  failure  times  graph  found  in  Dr.  Gularte’s 
article,  “Failure  of  Mechanically  Loaded  Laminated  Composites  Subjected  to 
Intense  Localized  Heating.” 

The  stresses  are  higher  near  the  boundaries  due  to  restrictions 
imposed  to  prevent  the  panel  from  expansion,  therefore  these  boundary 
conditions  or  restrictions  of  the  degrees  of  freedom  induce  higher  stress  away 
from  the  laser  beam  area. 

The  total  strain  energy  for  the  S.S.  case  was  326.7  pounds  per  inch. 

Laminated  composites  technology  to  improve  the  toughness  of 
laminated  composite  panels  is  limited.  Such  improvements  are  in  general 
usually  accompanied  by  a  sacrifice  in  other  properties  such  as  strength  at 
elevated  temperatures. 

Over  the  last  ten  years  or  so,  much  effort  was  concentrated  on- 
understanding  the  failure  mechanisms  and  predicting  the  strength  of 
compression  loaded  laminates.  The  fracture  of  composites  is  usually 
instantaneous  and  catastrophic,  which  makes  the  identification  of  critical 
failure  modes  hard  to  accomplish.  The  problem  is  further  complicated  because 
changes  in  material  properties,  or  the  presence  of  defects,  can  lead  to 
completely  different  failure  modes.  Hence,  an  analytical  model  accurate  for 
one  material  system  may  not  predict  failure  for  another  material  system.42 
Development  of  a  unified  model  which  can  be  applied  to  various  composite 

“Hahn  and  Williams,  Compression  Failure  Mechanics  in  Unidirectional  Composites,  p.  2. 
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material  systems  and  failure  modes  was  not  within  grasp.  However,  the 
current  approach  as  demonstrated  in  this  thesis,  can  identify  critical  failure 
modes  for  each  material  system  and  can  develop  a  unique  model  for  each 
failure  mode. 

Laminated  composites  performance  is  known  to  be  reduced  due  to 
the  substantial  amount  of  complex  load-induced  damage.  The  heterogeneity 
nature  of  the  laminated  composites  acts  not  just  as  a  crach  initiator,  but  also 
as  a  crack  arrester.43  Hence,  it  is  necessary  to  have  a  model  which  accounts 
for  the  average  effects  of  interests.  This  is  continuum  damage  mechanics, 
which  was  applied  successfully  to  isotropic  media  such  as  metal  and  concrete. 
However,  the  application  onto  laminated  orthotropic  composites  has  not  been 
widespread,  or  distinctly  successful. 

The  applications  objective  of  continuum  damage  mechanics  models 
to  laminated  composites  is  useful  for  engineering  design  purposes.  The 
alternative  would  be  to  attempt  to  solve  a  highly  anisotropic,  multiply 
connected,  non-linear,  non-symmetrical  boundary  value  problem. 

An  important  feature  that  stands  out  for  laminated  composites  is 
that  their  composite  tensile  modulus  is  greater  than  their  composite 
compressive  modulus.  Precise  theories  to  explain  this  feature  are  not 
available.  However,  one  possible  explanation  is  that  the  fibers  themselves  are 
less  stiff  in  compression  than  in  tension.  The  strain  hardening  observed  for 

“Allen,  Research  on  Damage  Models  for  Continuous  Fiber  Composites,  pp.  86-App.  7. 
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graphite/epoxy  composites  under  tension  is  attributed  to  improved  alignment 
of  internal  structure  of  the  graphite  fibers  during  loading.  Conversely,  initially 
imperfect  alignment  of  the  fiber  structure  may  grow  in  amplitude  during 
compression  loading.44 

As  mentioned  initially,  we  do  not  have  a  satisfactory  failure  criteria 
for  composite  materials  at  room  temperature  let  alone  elevated  temperatures. 

Composite  material  failure  criteria  will  have  to  account  for  the 
material  being  stronger  in  tension  than  in  compression.  The  elevated 
temperatures  effects  exaggerate  this,  i.e.  the  matrix  provides  lateral  resistant 
for  the  fibers,  and  as  the  matrix  is  heated  it  softens.  The  bottom  line  is  that 
with  current  fabrication  techniques,  irradiated  structures  are  markedly  more 
vulnerable  in  compression. 

The  results  of  this  work  have  been  studied  carefully  in  order  tc 
understand  the  compression  behavior  of  the  locally  heated  graphite/epoxy 
laminated  composite  panel.  The  predominant  failure  mode  has  been  identified 
as  shear  crippling,  due  to  the  high  bending  strains  in  the  fiber  in  the 
post-buckled  state.  This  leads  to  longitudinal  splitting  between  fibers  at  shear 
zone/s  because  of  the  required  kinematic  compatibility.  Shear  crippling  in 
laminated  composite  resembles  slip  lines  in  metals.46  The  most  important 
difficulty  in  predicting  failure  in  laminated  composite  panels  is  the  absence  of 


“Hahn  and  Williams,  Compression  Failure  Mechanics  in  Unidirectional  Composites,  p.  12. 
“Ibid,  pp.  16-17. 
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proven  failure  criteria  capable  of  handling  the  complex  loading  histories 
characteristic  of  these  responses.46  However,  the  buckling  response  in  this 
work  has  been  captured  reasonably  well. 

Equal  and  higher  laser  loading  rates,  as  that  used  in  this  work,  in 
conjunction  with  the  laminated  composite  panel  configuration,  will  almost 
always  result  in  catastrophic  failure. 

CONCLUSIONS: 

1.  Based  on  the  results  obtained,  the  analysis  program  of  the 
laser-induced  damage  calculation,  in  laminated  composites  panels,  of 
helicopters,  using  the  finite  element  method,  is  capable  of  predicting  the 
response  of  a  structure  under  compression  loading  subjected  to  laser  strike. 
Considering  the  mathematical  assumptions  made,  and  the  continuum  failure 
mechanism  approach  followed,  the  numerical  values  computed  should  be 
treated  as  suitable  for  the  conceptual  engineering  design  phase. 

2.  The  analysis  program  developed  in  this  thesis  met  the 
objectives  of  the  thesis  study,  i.e.  it  provided  enough  information  and  method 
for  the  prediction  of  helicopters  laminated  composites  panels  structural 
response  while  encountering  a  laser  strike.  The  analytical  simulations  have 

“Prantil,  Response  of  a  Thin  Cylindrical  Shell  Under  Lateral  Impulse  Loads,  pp.  57-59. 
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verified  the  computational  methods  for  studying  the  non-linear  buckling  failure 
response  of  compression  loaded,  localized  laser  heated,  laminated  composite 
panels. 


3.  Although  the  finite  element  simulations  used  here  represent  the 
state  of  the  art  capability  for  studying  laminated  composite  panels  response, 
however,  the  information  that  can  be  obtained  to  perform  an  accurate  laser 
lethality  loads  or  doses  assessment  is  limited. 

4.  It  is  concluded  from  the  stress  plots  that  at  higher  temperature 
(time  step  26),  although  stress  distribution  appears  to  have  the  same  pattern 
(as  indicated  by  the  stress  chart  bars  of  Figures  4  through  7  in  chapter  6),  the 
stress  quantities  were  considerably  higher. 

5.  With  current  fabrication  techniques  of  laminated  composites 
irradiated  structures  are  markedly  more  vulnerable  in  compression  than 
tension. 


SUGGESTIONS  FOR  FURTHER  WORKS: 


1.  Inclusion  of  spallation  effects  as  a  requirement  in  the  program 
output.  Incipient  spall  layer  can  be  one  model. 
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2.  Inclusion  of  ablation  effects  as  a  requirement  in  the  program 


output.  Different  variations  can  be  modeled. 

3.  Extension  of  the  present  work  to  models  that  account  for 
non-linear  mechanical  and  thermal  material  properties  in  each  of  the  three 
dimensions. 

4.  Inclusion  of  damping  and  investigating  its  effects. 

5.  inclusion  of  a  variety  of  laser  loading  and  beam  characteristics, 
together  with  different  target  geometries,  and  mechanical  and  thermal 
properties,  e.g.  pulse  versus  continuous  laser. 

6.  Investigation  of  the  compressive  shock  wave  damage. 

7.  The  magnitude  of  the  buckling  collapse  has  been  accurately 
predicted  here.  This  can  serve  as  a  relevant  input  to  any  future  suggested 
lethality  algorithm.  Reduction  in  buckling  threshold  can  represent  a  very 
viable  lethal  response  mechanism  to  be  incorporated  in  such  algorithm. 

8.  Physically  performing  accurate  detailed  experimentation 
program  under  laboratory  conditions  for  the  complete  and  variable 


- 107  - 


measurements  of  all  the  involved  parameters  in  the  laser-induced  damage  to 
mechanically  loaded  laminated  composite  panels  and  structures  used  in 
helicopters  and  aeroplanes. 
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